plus sc10x

load dependancies

load 10x data

GEX.seur <- readRDS("./GEX0811.seur.preAnno.rds")
GEX.seur
## An object of class Seurat 
## 18371 features across 25338 samples within 1 assay 
## Active assay: RNA (18371 features, 1229 variable features)
##  3 dimensional reductions calculated: pca, tsne, umap
color.FB <- ggsci::pal_igv("default")(49)[c(2,13,33,1,15,
                                            34,26,28,40,18)]

color.cnt <- color.FB[c(10,9,7,
                        5,4,1)]
DimPlot(GEX.seur, reduction = "umap", label = T, group.by = "preAnno") +
  DimPlot(GEX.seur, reduction = "umap", label = F, group.by = "cnt", cols = color.cnt)

DimPlot(GEX.seur, reduction = "umap", label = T,label.size = 2.4, group.by = "preAnno1") +
  DimPlot(GEX.seur, reduction = "umap", label = F, group.by = "preAnno2")

VlnPlot(GEX.seur, features = c("nFeature_RNA", "nCount_RNA", "percent.mt","percent.rb"), ncol = 2, pt.size = 0.1, group.by = "preAnno1")

VlnPlot(GEX.seur, features = c("nFeature_RNA", "nCount_RNA", "percent.mt","percent.rb"), ncol = 2, pt.size = 0.1, group.by = "preAnno1", split.by = "DoubletFinder0.05")
## The default behaviour of split.by has changed.
## Separate violin plots are now plotted side-by-side.
## To restore the old behaviour of a single split violin,
## set split.plot = TRUE.
##       
## This message will be shown once per session.

VlnPlot(GEX.seur, features = c("nFeature_RNA", "nCount_RNA", "percent.mt","percent.rb"), ncol = 2, pt.size = 0.1, group.by = "FB.info")

VlnPlot(GEX.seur, features = c("nFeature_RNA", "nCount_RNA", "percent.mt","percent.rb"), ncol = 2, pt.size = 0.1, group.by = "FB.info", split.by = "DoubletFinder0.05")

advanced filtering

only keep DF0.05 Singlets
remove C1

GEX.seur <- subset(GEX.seur, subset = DoubletFinder0.05=="Singlet" & seurat_clusters %in% setdiff(c(0:19),c(15,18,19)))
GEX.seur <- subset(GEX.seur, subset = nFeature_RNA < 3000 & nCount_RNA < 8000)
GEX.seur
## An object of class Seurat 
## 18371 features across 23692 samples within 1 assay 
## Active assay: RNA (18371 features, 1229 variable features)
##  3 dimensional reductions calculated: pca, tsne, umap
VlnPlot(GEX.seur, features = c("nFeature_RNA", "nCount_RNA", "percent.mt","percent.rb"), ncol = 2, pt.size = 0.1, group.by = "preAnno1")

VlnPlot(GEX.seur, features = c("nFeature_RNA", "nCount_RNA", "percent.mt","percent.rb"), ncol = 2, pt.size = 0, group.by = "preAnno1")

VlnPlot(GEX.seur, features = c("nFeature_RNA", "nCount_RNA", "percent.mt","percent.rb"), ncol = 2, pt.size = 0.1, group.by = "FB.info")

VlnPlot(GEX.seur, features = c("nFeature_RNA", "nCount_RNA", "percent.mt","percent.rb"), ncol = 2, pt.size = 0, group.by = "FB.info")

plota <- FeatureScatter(GEX.seur, feature1 = "nCount_RNA", feature2 = "percent.mt", group.by = "cnt", cols = color.cnt) 
plotb <- FeatureScatter(GEX.seur, feature1 = "nCount_RNA", feature2 = "nFeature_RNA", group.by = "cnt", cols = color.cnt) 
plotc <- FeatureScatter(GEX.seur, feature1 = "nCount_RNA", feature2 = "percent.rb", group.by = "cnt", cols = color.cnt) 
plota + plotb + plotc

GEX.seur$FB.new <- factor(as.character(GEX.seur$FB.info),
                          levels = setdiff(levels(GEX.seur$FB.info),
                                           c("Doublet","Negative"))[c(10,8:9,6:7,
                                                                      5,3:4,1:2)])
color.FBnew <- color.FB[c(10,8:9,6:7,5,3:4,1:2)]
par(mar=c(6,4,2,3))
sl_stat <- table(GEX.seur$FB.new)
barplot(sl_stat,ylim = c(0,4200),
        #col = c("#FF6C91","lightgrey",color.FB),
        col = c(color.FBnew),
        main = "Feature Barcode statistics",cex.names = 0.75, xaxt = "n")
axis(1,1:10*1.2-0.48 ,levels(GEX.seur@meta.data$FB.new), las=3, cex.axis=0.85)
text(x=1:10*1.2-0.45,y=sl_stat+285,paste0(sl_stat,"\n",100*round(as.numeric(sl_stat/sum(sl_stat)),4),"%"),cex = 0.75)

new

re-clustering

GEX.seur <- FindVariableFeatures(GEX.seur, selection.method = "vst", nfeatures = 1800)

# Identify the 10 most highly variable genes
#top10 <- head(VariableFeatures(GEX.seur), 10)
top20 <- head(VariableFeatures(GEX.seur), 20)

# plot variable features with and without labels
plot1 <- VariableFeaturePlot(GEX.seur)
plot2 <- LabelPoints(plot = plot1, points = top20, repel = T)
plot1 + plot2

head(VariableFeatures(GEX.seur), 300)
##   [1] "Spp1"          "Ccl4"          "Cxcl10"        "Ccl5"         
##   [5] "Ccl3"          "Spp1-EGFP"     "Ccl12"         "Hist1h1b"     
##   [9] "Mki67"         "Cxcl9"         "Ifi27l2a"      "Ccl2"         
##  [13] "Cxcl13"        "Apoe"          "Rsad2"         "Top2a"        
##  [17] "Xist"          "Fn1"           "Cst7"          "Hist1h2ae"    
##  [21] "Cenpf"         "Lpl"           "Ifitm3"        "Hmgb2"        
##  [25] "Fabp5"         "Hist1h2ap"     "Lgals3"        "Ifit3"        
##  [29] "Ch25h"         "Ifit2"         "Stmn1"         "Prc1"         
##  [33] "Iigp1"         "Isg15"         "Ifit1"         "Pclaf"        
##  [37] "Ccl7"          "Cd52"          "Hmmr"          "Ube2c"        
##  [41] "Egr1"          "Gpnmb"         "Nusap1"        "Mmp12"        
##  [45] "Cenpa"         "Cd72"          "Cd69"          "Il1b"         
##  [49] "Apoc1"         "Csf1"          "Gm26737"       "Aspm"         
##  [53] "Cd9"           "Ctsd"          "Cenpe"         "Birc5"        
##  [57] "Cd74"          "H2-K1"         "Cdca8"         "Pbld1"        
##  [61] "Ifi204"        "Gdf15"         "Lgals3bp"      "H2-D1"        
##  [65] "H2afz"         "Igf1"          "Lyz2"          "Hist1h3c"     
##  [69] "Tpx2"          "Ccl6"          "Slpi"          "Ctsb"         
##  [73] "Saa3"          "B230312C02Rik" "Arg1"          "Il12b"        
##  [77] "Fth1"          "Bst2"          "Serpina3g"     "Smc2"         
##  [81] "Mif"           "Kif11"         "Lgals1"        "Acod1"        
##  [85] "Ccl8"          "Esco2"         "Fxyd5"         "Ttr"          
##  [89] "B2m"           "Smc4"          "H2afx"         "Knl1"         
##  [93] "Hist1h1e"      "Ftl1"          "Ccnb1"         "Cd5l"         
##  [97] "Cdkn1a"        "Ms4a7"         "Cd63"          "Pbk"          
## [101] "Rrm2"          "Mrc1"          "Hmox1"         "Igfbp5"       
## [105] "Ifi211"        "Slfn5"         "Cybb"          "Hba-a2"       
## [109] "Hba-a1"        "Kif15"         "Gm26885"       "Dqx1"         
## [113] "Rcan1"         "Hist2h2ac"     "Wfdc17"        "Kif23"        
## [117] "Ctsz"          "Mis18bp1"      "Lilrb4a"       "Ifi207"       
## [121] "Cdca3"         "Kif20b"        "Ccnb2"         "Apoc2"        
## [125] "Mt1"           "4933407L21Rik" "3110039M20Rik" "Gbp2"         
## [129] "Syngr1"        "H2-Q7"         "Atad2"         "Ifi203"       
## [133] "Ctsl"          "Ccr1"          "Enpp2"         "Marco"        
## [137] "Gm15987"       "Pf4"           "Gm42047"       "Cmpk2"        
## [141] "Cdc20"         "Ifi44"         "Ifit3b"        "Bub1b"        
## [145] "Ccl9"          "Cxcl2"         "Tnfsf18"       "Anln"         
## [149] "Racgap1"       "Dennd5b"       "Aldh1a3"       "Gfod2"        
## [153] "Wdcp"          "Atf3"          "Sgo2a"         "Tspo"         
## [157] "Ifitm1"        "Ckap2l"        "Tc2n"          "Cspg4"        
## [161] "Spc24"         "Ndc80"         "Gm43409"       "Plac8"        
## [165] "Irf7"          "Tubb4b"        "Oasl2"         "Ctla2a"       
## [169] "Rbms3"         "Myh4"          "Lmnb1"         "C4b"          
## [173] "Gnas"          "H2-Ab1"        "Hist1h4d"      "2210409D07Rik"
## [177] "Sh3tc2"        "Cdk1"          "Cox6a2"        "Gapdh"        
## [181] "Gm12863"       "Dcx"           "Timp2"         "Kcnq1ot1"     
## [185] "Phf11b"        "Spc25"         "Ly6a"          "Ccl24"        
## [189] "Ank"           "Hist1h2ab"     "Cdca2"         "Ly86"         
## [193] "Baiap2l2"      "Fos"           "Tk1"           "Tsix"         
## [197] "Lig1"          "Cks2"          "Cstb"          "Ifi213"       
## [201] "Clspn"         "Ccnf"          "Nfasc"         "Gfap"         
## [205] "Lockd"         "Cks1b"         "Cp"            "Snca"         
## [209] "Gm4951"        "Cxcl14"        "Tex11"         "Fam20c"       
## [213] "Ckap2"         "Parp14"        "Cxcl16"        "Lrrc31"       
## [217] "Gm10457"       "Tmpo"          "Ccna2"         "Ifi209"       
## [221] "Plk1"          "Rbpms2"        "Lsamp"         "Ifi205"       
## [225] "Loxl2"         "Vim"           "Tnr"           "Prokr1"       
## [229] "Ifi208"        "Apob"          "Ifnb1"         "Fbxo5"        
## [233] "Ncapg2"        "Tlr2"          "Incenp"        "Pcna"         
## [237] "Gm8251"        "mt-Nd1"        "Diaph3"        "Ncapg"        
## [241] "Zbp1"          "Capg"          "Emp3"          "Ndnf"         
## [245] "Gm35558"       "Rnd3"          "Rad51ap1"      "F830016B08Rik"
## [249] "Apoc4"         "Mamdc2"        "Tubb5"         "Aurkb"        
## [253] "Trim59"        "Crybb1"        "Gadd45b"       "Srgn"         
## [257] "Clec4d"        "Olfr433"       "Fgl2"          "Egr3"         
## [261] "Gas2l3"        "Kif4"          "Nr4a3"         "Il1rn"        
## [265] "Slfn1"         "Sgo1"          "Stat1"         "Gm43568"      
## [269] "Cd83"          "Oasl1"         "Ifi206"        "Gpr84"        
## [273] "Tuba1b"        "Usp18"         "Id2"           "Dtl"          
## [277] "Dlgap5"        "Mt2"           "Arl6ip1"       "Anxa5"        
## [281] "H2-Aa"         "mt-Co1"        "Hist1h3e"      "Ccdc34"       
## [285] "Rgs17"         "Olfr1369-ps1"  "Rnase2a"       "Cpd"          
## [289] "E2f8"          "Syt1"          "Chp2"          "Fam57b"       
## [293] "Ctnna3"        "Cbx5"          "Rrm1"          "Gm29773"      
## [297] "Gdf3"          "Nuf2"          "Herc6"         "Ier3"
# exclude MT genes  and more 
#   add sex-related Xist/Tsix
DIG <- grep("^Tra|^Trb|^Trg|^Trd|^Tcr|^Igm|^Igh|^Igk|^Igl|Jchain|Mzb1|Vpreb|Lars2|Jun|Fos|Egr1|^Hsp|^Rps|^Rpl|Hbb-|Hba-|^Dnaj|^AC|^AI|^AA|^AW|^AY|^BC|^Gm|^Hist|Rik$|-ps|Xist|Tsix|^Ifi|^Isg|^Mcm",
            rownames(GEX.seur),value = T)
CC_gene <- Hmisc::capitalize(tolower(as.vector(unlist(cc.genes.updated.2019))))
MT_gene <- grep("^mt-",rownames(GEX.seur),value = T)

VariableFeatures(GEX.seur) <- setdiff(VariableFeatures(object = GEX.seur),
                                                c(MT_gene,
                                                  DIG, 
                                                  CC_gene) )

GEX.seur <- RunPCA(GEX.seur, features = VariableFeatures(GEX.seur), seed.use = 868)
## PC_ 1 
## Positive:  Cst7, Cd72, B2m, Ch25h, Lgals3bp, Cd52, Fth1, Ctsb, Ctsz, Ccl3 
##     Cd9, Bst2, Fxyd5, Ctsd, H2-K1, Mif, Lpl, H2-D1, Fabp5, Spp1-EGFP 
##     Lilrb4a, Csf1, Tspo, Ccl2, Ccl4, Cd63, Lyz2, Zbp1, Oasl2, Fn1 
## Negative:  Crybb1, Ccr5, Sox4, Fchsd2, Hpgd, Clec4a3, Bank1, Fcgrt, Lst1, Klf7 
##     Mrc1, Igfbp4, Serpinf1, Camk2n1, Rgs7bp, Fry, Gpr165, Cbfa2t3, Tnfrsf17, Il7r 
##     Cryba4, Gbp7, Filip1l, Lrba, Ppm1l, Rdx, Zbtb20, Ccr1, Ddit4, Nav2 
## PC_ 2 
## Positive:  Prc1, Pclaf, Pbk, Knl1, Esco2, Kif15, Aspm, Spc24, Stmn1, Neil3 
##     Smc2, Racgap1, Mis18bp1, Lockd, Bub1b, Ccna2, Sgo2a, Ccnb1, Ccnf, H2afx 
##     Plk1, Tk1, Cit, Spc25, Ncapg, Diaph3, Shcbp1, Trim59, Kif4, Fbxo5 
## Negative:  Ctsd, Ctsb, Tyrobp, B2m, H2-K1, Cst7, Cd9, Timp2, Fth1, H2-D1 
##     Cd63, Ctsz, Apoe, Lgals3bp, Lyz2, Ctsl, Ftl1, Grn, Npc2, Mpeg1 
##     Slfn5, Ly86, Ccl3, H2-Q7, H2-T23, Bst2, Oasl2, Syngr1, Cd72, Ccl4 
## PC_ 3 
## Positive:  Rsad2, Iigp1, Slfn5, Oasl2, Parp14, Stat1, Usp18, Fgl2, Phf11d, Cxcl10 
##     Rtp4, Herc6, Slfn8, Trim30a, Ly6e, Ly6a, Slfn2, 9930111J21Rik2, Sp100, Oasl1 
##     Irf7, Rnf213, Zup1, Gbp2, Stat2, Tor3a, Xaf1, Ddx58, Ccl12, Etnk1 
## Negative:  Tmsb4x, Ftl1, Lpl, Tyrobp, Igf1, Mt1, Fabp5, Nme2, Serf2, Mif 
##     Rbm3, Myl6, Hint1, Cstb, Nme1, Spp1-EGFP, Gapdh, Lgals3, Cox6c, Atox1 
##     Emp3, Gpx1, Atp5g1, Gng5, Pkm, Gas6, Scd2, Calr, Fxyd5, Aldoa 
## PC_ 4 
## Positive:  Ctsd, Cd9, Pmp22, Nav2, Samsn1, Fam102b, Hpgd, Csmd3, Lrba, Sgk1 
##     Ctsz, Cadm1, Tent5c, Syngr1, Thrsp, Pdgfb, Srgn, Lag3, Ccl6, Gpr183 
##     Cd34, Il7r, Ccl9, Prc1, Cd14, Ccnb1, Aspm, Timp2, Fth1, Nav3 
## Negative:  Tmsb4x, Mrc1, Rbm3, Rsad2, Iigp1, Ccl12, Igfbp4, Slfn5, Crybb1, Cxcl10 
##     Fcgrt, Gas6, Gbp7, Nme2, Cryba4, Oasl1, Oasl2, Usp18, Chd4, Gbp2 
##     Eif3a, Parp14, Ptma, Ccnd1, Lst1, Eif2ak2, Tpr, Hint1, Zbp1, Zfas1 
## PC_ 5 
## Positive:  Tmsb4x, Apoe, Ccnb1, Aspm, Lyz2, Knstrn, Tyrobp, Cdkn3, Ftl1, Plk1 
##     Cep55, Ctsl, Zfas1, Cd52, Racgap1, Cox8a, Fth1, Apoc1, Lockd, Gpx1 
##     Prc1, H2afx, Bub1b, Tuba1c, B2m, Selenoh, Ramp1, Kif20a, Cox6c, Hacd4 
## Negative:  Lig1, Dhfr, Dut, Macf1, Atad5, Fignl1, Ccl4, Cdt1, Fam102b, Nsd2 
##     Dnmt1, Nav3, Wdhd1, Rif1, Kcnq1ot1, Rad54l, Rmi2, Topbp1, Timeless, Chaf1a 
##     Brca2, Ncapg2, Pmp22, Ccnd1, Tcf19, Pdgfb, Brca1, Nop56, Hnrnpab, Lilrb4a
length(VariableFeatures(GEX.seur))
## [1] 1389
head(VariableFeatures(GEX.seur),300)
##   [1] "Spp1"      "Ccl4"      "Cxcl10"    "Ccl5"      "Ccl3"      "Spp1-EGFP"
##   [7] "Ccl12"     "Cxcl9"     "Ccl2"      "Cxcl13"    "Apoe"      "Rsad2"    
##  [13] "Fn1"       "Cst7"      "Lpl"       "Fabp5"     "Lgals3"    "Ch25h"    
##  [19] "Stmn1"     "Prc1"      "Iigp1"     "Pclaf"     "Ccl7"      "Cd52"     
##  [25] "Gpnmb"     "Mmp12"     "Cd72"      "Cd69"      "Il1b"      "Apoc1"    
##  [31] "Csf1"      "Aspm"      "Cd9"       "Ctsd"      "Cd74"      "H2-K1"    
##  [37] "Pbld1"     "Gdf15"     "Lgals3bp"  "H2-D1"     "H2afz"     "Igf1"     
##  [43] "Lyz2"      "Ccl6"      "Slpi"      "Ctsb"      "Saa3"      "Arg1"     
##  [49] "Il12b"     "Fth1"      "Bst2"      "Serpina3g" "Smc2"      "Mif"      
##  [55] "Lgals1"    "Acod1"     "Ccl8"      "Esco2"     "Fxyd5"     "Ttr"      
##  [61] "B2m"       "H2afx"     "Knl1"      "Ftl1"      "Ccnb1"     "Cd5l"     
##  [67] "Cdkn1a"    "Ms4a7"     "Cd63"      "Pbk"       "Mrc1"      "Hmox1"    
##  [73] "Igfbp5"    "Slfn5"     "Cybb"      "Kif15"     "Dqx1"      "Rcan1"    
##  [79] "Wfdc17"    "Ctsz"      "Mis18bp1"  "Lilrb4a"   "Apoc2"     "Mt1"      
##  [85] "Gbp2"      "Syngr1"    "H2-Q7"     "Ctsl"      "Ccr1"      "Enpp2"    
##  [91] "Marco"     "Pf4"       "Cmpk2"     "Bub1b"     "Ccl9"      "Cxcl2"    
##  [97] "Tnfsf18"   "Racgap1"   "Dennd5b"   "Aldh1a3"   "Gfod2"     "Wdcp"     
## [103] "Atf3"      "Sgo2a"     "Tspo"      "Tc2n"      "Cspg4"     "Spc24"    
## [109] "Plac8"     "Irf7"      "Oasl2"     "Ctla2a"    "Rbms3"     "Myh4"     
## [115] "Lmnb1"     "C4b"       "Gnas"      "H2-Ab1"    "Sh3tc2"    "Cox6a2"   
## [121] "Gapdh"     "Dcx"       "Timp2"     "Kcnq1ot1"  "Phf11b"    "Spc25"    
## [127] "Ly6a"      "Ccl24"     "Ank"       "Ly86"      "Baiap2l2"  "Tk1"      
## [133] "Lig1"      "Cstb"      "Ccnf"      "Nfasc"     "Gfap"      "Lockd"    
## [139] "Cp"        "Snca"      "Cxcl14"    "Tex11"     "Fam20c"    "Parp14"   
## [145] "Cxcl16"    "Lrrc31"    "Ccna2"     "Plk1"      "Rbpms2"    "Lsamp"    
## [151] "Loxl2"     "Vim"       "Tnr"       "Prokr1"    "Apob"      "Ifnb1"    
## [157] "Fbxo5"     "Ncapg2"    "Tlr2"      "Incenp"    "Diaph3"    "Ncapg"    
## [163] "Zbp1"      "Capg"      "Emp3"      "Ndnf"      "Rnd3"      "Apoc4"    
## [169] "Mamdc2"    "Tubb5"     "Trim59"    "Crybb1"    "Gadd45b"   "Srgn"     
## [175] "Clec4d"    "Olfr433"   "Fgl2"      "Egr3"      "Kif4"      "Nr4a3"    
## [181] "Il1rn"     "Slfn1"     "Sgo1"      "Stat1"     "Cd83"      "Oasl1"    
## [187] "Gpr84"     "Tuba1b"    "Usp18"     "Id2"       "Mt2"       "Arl6ip1"  
## [193] "Anxa5"     "H2-Aa"     "Ccdc34"    "Rgs17"     "Rnase2a"   "Cpd"      
## [199] "Syt1"      "Chp2"      "Fam57b"    "Ctnna3"    "Gdf3"      "Herc6"    
## [205] "Ier3"      "Phf11d"    "Cdk15"     "Zdbf2"     "Otor"      "Cldn11"   
## [211] "Sox2ot"    "Lpar1"     "Akr1b7"    "Tril"      "Fancd2os"  "Fgf6"     
## [217] "Vmn2r30"   "Osbpl5"    "Fat2"      "Asic2"     "Cyb561"    "Sox9"     
## [223] "F2rl1"     "Blk"       "Igsf11"    "Ripply3"   "Syt4"      "Crtac1"   
## [229] "Col4a5"    "Cdkn3"     "Lhfpl4"    "Hif1a"     "Knstrn"    "Fabp3"    
## [235] "H2afv"     "Clmn"      "Cst6"      "Cck"       "Ldha"      "Neil3"    
## [241] "Rnf213"    "Slfn2"     "Selenow"   "Itgax"     "Tpi1"      "Ptger4"   
## [247] "Alas2"     "Il4i1"     "Fcgr4"     "Stil"      "Ezh2"      "Ms4a4c"   
## [253] "Klhl33"    "Tuba1c"    "Pkm"       "Hirip3"    "Ecrg4"     "Igfbp2"   
## [259] "Il12rb1"   "Cd24a"     "Cd40"      "Rgs1"      "Olfr889"   "Dut"      
## [265] "Rab7b"     "Rtp4"      "Kmo"       "Mcf2l"     "Olr1"      "Meg3"     
## [271] "E2f7"      "Anp32b"    "Folr2"     "Troap"     "Adam33"    "Cit"      
## [277] "Gbp5"      "Ly6e"      "Shcbp1"    "Melk"      "Plau"      "Zfp9"     
## [283] "Oas3"      "Colec12"   "Cdkn2c"    "Sag"       "Otos"      "Adcy10"   
## [289] "Ermn"      "Pde1a"     "Grin2b"    "Hs3st4"    "Tex15"     "Rasd2"    
## [295] "Tox3"      "Ccdc151"   "Ptprk"     "Ppp1r27"   "Dpf3"      "Sycp2l"
DimPlot(GEX.seur, reduction = "pca",dims = 1:2, group.by = "orig.ident") +
  DimPlot(GEX.seur, reduction = "pca",dims = 3:4, group.by = "orig.ident")

DimPlot(GEX.seur, reduction = "pca",dims = 1:2, group.by = "FB.new", cols = color.FBnew, ) +
  DimPlot(GEX.seur, reduction = "pca",dims = 3:4, group.by = "FB.new", cols = color.FBnew)

DimPlot(GEX.seur, reduction = "pca",dims = 1:2, group.by = "cnt", cols = color.cnt) +
  DimPlot(GEX.seur, reduction = "pca",dims = 3:4, group.by = "cnt", cols = color.cnt)

cowplot::plot_grid(
FeaturePlot(GEX.seur, features = "Spp1", reduction = "pca",dims = 1:2),
  DimPlot(GEX.seur, reduction = "pca",dims = 1:2, group.by = "cnt", cols = color.cnt),
rel_widths = c(5,5.6),
ncol = 2)

cowplot::plot_grid(
FeaturePlot(GEX.seur, features = "Spp1-EGFP", reduction = "pca",dims = 1:2),
  DimPlot(GEX.seur, reduction = "pca",dims = 1:2, group.by = "cnt", cols = color.cnt),
rel_widths = c(5,5.6),
ncol = 2)

DimPlot(GEX.seur, reduction = "pca",dims = 1:2, group.by = "FB.new", split.by = "FB.new", cols = color.FBnew, ncol = 5)

DimPlot(GEX.seur, reduction = "pca",dims = 1:2, group.by = "cnt", split.by = "cnt", cols = color.cnt, ncol = 3)

FeaturePlot(GEX.seur, features = "Spp1", reduction = "pca",dims = 1:2, split.by = "cnt")

FeaturePlot(GEX.seur, features = "Spp1-EGFP", reduction = "pca",dims = 1:2, split.by = "cnt")

FeaturePlot(GEX.seur, features = "P2ry12", reduction = "pca",dims = 1:2, split.by = "cnt")

FeaturePlot(GEX.seur, features = "Mki67", reduction = "pca",dims = 1:2, split.by = "cnt")

cowplot::plot_grid(
FeaturePlot(subset(GEX.seur,subset=cnt %in% c("P05.CTR","P05.MIG","P05.GFP")), features = "Spp1", reduction = "pca",dims = 1:2, split.by = "cnt", ) ,#& xlim(c(-7.5,33)) & ylim(c(-11,49)) ,
  FeaturePlot(subset(GEX.seur,subset=cnt %in% c("P28.CTR","P28.MIG","P28.GFP")), features = "Spp1", reduction = "pca",dims = 1:2, split.by = "cnt") ,#& xlim(c(-7.5,33)) & ylim(c(-11,49)),
ncol = 1)

cowplot::plot_grid(
FeaturePlot(subset(GEX.seur,subset=cnt %in% c("P05.CTR","P05.MIG","P05.GFP")), features = "Spp1-EGFP", reduction = "pca",dims = 1:2, split.by = "cnt") ,
  FeaturePlot(subset(GEX.seur,subset=cnt %in% c("P28.CTR","P28.MIG","P28.GFP")), features = "Spp1-EGFP", reduction = "pca",dims = 1:2, split.by = "cnt") ,
ncol = 1)

cowplot::plot_grid(
FeaturePlot(GEX.seur, features = "P2ry12", reduction = "pca",dims = 1:2),
  DimPlot(GEX.seur, reduction = "pca",dims = 1:2, group.by = "cnt", cols = color.cnt),
rel_widths = c(5,5.6),
ncol = 2)

cowplot::plot_grid(
FeaturePlot(subset(GEX.seur,subset=cnt %in% c("P05.CTR","P05.MIG","P05.GFP")), features = "P2ry12", reduction = "pca",dims = 1:2, split.by = "cnt") ,
  FeaturePlot(subset(GEX.seur,subset=cnt %in% c("P28.CTR","P28.MIG","P28.GFP")), features = "P2ry12", reduction = "pca",dims = 1:2, split.by = "cnt") ,
ncol = 1)

cowplot::plot_grid(
FeaturePlot(GEX.seur, features = "Mki67", reduction = "pca",dims = 1:2),
  DimPlot(GEX.seur, reduction = "pca",dims = 1:2, group.by = "cnt", cols = color.cnt),
rel_widths = c(5,5.6),
ncol = 2)

cowplot::plot_grid(
FeaturePlot(subset(GEX.seur,subset=cnt %in% c("P05.CTR","P05.MIG","P05.GFP")), features = "Mki67", reduction = "pca",dims = 1:2, split.by = "cnt") ,
  FeaturePlot(subset(GEX.seur,subset=cnt %in% c("P28.CTR","P28.MIG","P28.GFP")), features = "Mki67", reduction = "pca",dims = 1:2, split.by = "cnt") ,
ncol = 1)

DimHeatmap(GEX.seur, dims = 1:12, cells = 6000, balanced = TRUE,ncol = 4,nfeatures = 80)

(GEX.seur@reductions$pca@feature.loadings %>% as.data.frame)[,1:8] %>% arrange(desc(PC_1) ) %>% head(40)
##                 PC_1        PC_2         PC_3          PC_4         PC_5
## Cst7      0.13033723 -0.04914273 -0.024344645  0.0391614600  0.022557581
## Cd72      0.10462524 -0.03082428 -0.025466799 -0.0060531634 -0.020052539
## B2m       0.10373478 -0.05260759  0.042947295  0.0302277177  0.060053627
## Ch25h     0.10215044 -0.02902052 -0.016884600 -0.0136793196 -0.058608146
## Lgals3bp  0.10039744 -0.04596373  0.067935444  0.0083731503  0.024409626
## Cd52      0.10030771 -0.03034810 -0.012582220 -0.0374747502  0.072467231
## Fth1      0.10017184 -0.04848443 -0.014527262  0.0583658027  0.070162498
## Ctsb      0.09965710 -0.05693450 -0.041613614  0.0405126748  0.037755196
## Ctsz      0.09570110 -0.04719809 -0.016882643  0.0814963494  0.019413379
## Ccl3      0.09494543 -0.03390996 -0.033316477  0.0319485657 -0.066979156
## Cd9       0.09446644 -0.04879779 -0.007663390  0.1393953163 -0.020040494
## Bst2      0.09385601 -0.03219153  0.043282808 -0.0558837565  0.045345922
## Fxyd5     0.09371542 -0.01947810 -0.057627566  0.0095793204  0.004101436
## Ctsd      0.09312216 -0.07980194  0.010063702  0.1698541195 -0.012359922
## H2-K1     0.09307069 -0.04946374  0.061854974  0.0557680327  0.040462048
## Mif       0.09164890 -0.01956187 -0.072324527  0.0098647201 -0.003255225
## Lpl       0.08980026 -0.01488866 -0.090240695 -0.0003338312 -0.040050573
## H2-D1     0.08873270 -0.04841797  0.060894356  0.0517284292  0.055801786
## Fabp5     0.08762805 -0.02600469 -0.079178279  0.0150022563 -0.010474565
## Spp1-EGFP 0.08739725 -0.02331045 -0.066788887  0.0056088220 -0.017086488
## Lilrb4a   0.08727945 -0.02537670 -0.048874340  0.0171707711 -0.067661249
## Csf1      0.08626413 -0.02907911 -0.024130726  0.0133250495 -0.058306480
## Tspo      0.08615673 -0.02259411  0.020394771 -0.0420595092  0.055264464
## Ccl2      0.08516958 -0.02556872  0.044195824 -0.0586610844 -0.033474700
## Ccl4      0.08508932 -0.03059175 -0.024160255  0.0181325362 -0.085269328
## Cd63      0.08310536 -0.04724694 -0.054186027  0.0300100818  0.030584581
## Lyz2      0.08260736 -0.03988305 -0.022641236  0.0189047277  0.089614900
## Zbp1      0.08228250 -0.02797594  0.070085144 -0.0694688088  0.013992002
## Oasl2     0.08192573 -0.03214081  0.116312337 -0.0799565678  0.009603802
## Fn1       0.08104928 -0.02006132 -0.040377111 -0.0127196754 -0.013037203
## Irf7      0.08059542 -0.02868662  0.086758150 -0.0685628905  0.015186291
## Apoe      0.07978716 -0.04636109 -0.025883130 -0.0538887914  0.137889836
## H2-Q7     0.07970411 -0.03245787  0.052403842  0.0215995009  0.030129185
## Rsad2     0.07945714 -0.02853418  0.127246449 -0.1123588025 -0.011752638
## Spp1      0.07916311 -0.02501669 -0.056085597  0.0140606130 -0.004654707
## Cxcl16    0.07808397 -0.02833363 -0.001287169  0.0348527122 -0.019926816
## Phf11b    0.07709413 -0.02253058  0.070077143 -0.0586345152  0.009170246
## Ftl1      0.07580178 -0.03755221 -0.113949095 -0.0174402270  0.084352113
## Gapdh     0.07571862  0.00168726 -0.066102735 -0.0014192182 -0.005445161
## Slfn5     0.07523311 -0.03439980  0.121948327 -0.0986415897  0.015533179
##                   PC_6         PC_7         PC_8
## Cst7      -0.044944616  0.006920554  0.037016038
## Cd72       0.025995912 -0.028995407  0.075412988
## B2m       -0.117409096  0.028275984 -0.024816009
## Ch25h      0.091269216 -0.029847492  0.101975365
## Lgals3bp  -0.090008619  0.054464751 -0.040151864
## Cd52      -0.070996534 -0.009698901  0.013347348
## Fth1      -0.088306196 -0.045837095  0.002304896
## Ctsb      -0.057557479  0.145524810 -0.020520516
## Ctsz      -0.059787151  0.004880728 -0.053875683
## Ccl3       0.092608183  0.013953343  0.110391444
## Cd9       -0.034012305 -0.023169678 -0.038052392
## Bst2      -0.044777599 -0.015277668 -0.036204519
## Fxyd5      0.020256812 -0.047032567  0.025598323
## Ctsd      -0.073937168  0.024783903 -0.049834277
## H2-K1     -0.098054435  0.023869165 -0.004797431
## Mif       -0.006587504 -0.031618471 -0.033623015
## Lpl        0.060525346  0.027615059  0.042608584
## H2-D1     -0.119351346  0.042041655 -0.001854407
## Fabp5      0.028043129  0.026398805  0.049128937
## Spp1-EGFP  0.070929233 -0.050363638  0.079791275
## Lilrb4a    0.088976630  0.015712180  0.104512794
## Csf1       0.076990407  0.029524221  0.057088988
## Tspo      -0.038027118 -0.051554132 -0.019477712
## Ccl2       0.085629727 -0.054186904  0.111073502
## Ccl4       0.127044371 -0.011092153  0.146989600
## Cd63      -0.057374763  0.093133316 -0.008161522
## Lyz2      -0.120994907  0.120351793 -0.018108488
## Zbp1       0.020027446 -0.046656939  0.018233516
## Oasl2      0.001815700 -0.003693521 -0.022303578
## Fn1        0.073433024 -0.049732703  0.067528229
## Irf7      -0.001806470 -0.035682258  0.009648752
## Apoe      -0.110169908  0.213888729  0.111025894
## H2-Q7     -0.069706644  0.005301151  0.009538434
## Rsad2      0.078424509 -0.042403429  0.053046329
## Spp1       0.054218209 -0.018373564  0.085418704
## Cxcl16     0.006855720 -0.012366960  0.035189155
## Phf11b     0.025043245 -0.072215760  0.010316775
## Ftl1      -0.036603625  0.027822457 -0.030990694
## Gapdh     -0.011628249 -0.042104589 -0.070552164
## Slfn5      0.025909198  0.014281558 -0.039929855
(GEX.seur@reductions$pca@feature.loadings %>% as.data.frame)[,1:8] %>% arrange(PC_1 ) %>% head(40)
##                 PC_1          PC_2          PC_3         PC_4         PC_5
## Crybb1   -0.06463817  0.0294974727 -0.0354208164 -0.097408764  0.006862056
## Ccr5     -0.06204224  0.0110291106  0.0140474164 -0.043660322 -0.028842188
## Sox4     -0.04764682  0.0121422496 -0.0171804830 -0.006418549 -0.041590884
## Fchsd2   -0.04629952  0.0227555136 -0.0045365074 -0.049839575 -0.046466318
## Hpgd     -0.03950759 -0.0059255018  0.0374309333  0.088391440 -0.037506101
## Clec4a3  -0.03504186 -0.0010099934  0.0138714112 -0.002062804  0.005103966
## Bank1    -0.03184373  0.0124236882 -0.0018493449 -0.016879860 -0.037725693
## Fcgrt    -0.03149697  0.0198491582 -0.0373920652 -0.090065353  0.011228072
## Lst1     -0.03141994  0.0077358960 -0.0143684268 -0.070862561  0.049164853
## Klf7     -0.03067361  0.0135257120 -0.0053927050 -0.030175833 -0.019121449
## Mrc1     -0.03063504  0.0457878200 -0.0460962423 -0.119637080  0.036337511
## Igfbp4   -0.02825798  0.0380816918 -0.0362660542 -0.104365500  0.021842448
## Serpinf1 -0.02810006  0.0142248546 -0.0151006007 -0.037692334  0.009542200
## Camk2n1  -0.02488802  0.0103614022  0.0006635641 -0.010470197 -0.026282061
## Rgs7bp   -0.02397541  0.0188955579 -0.0266946601 -0.039995923 -0.006532197
## Fry      -0.02292624  0.0080681519  0.0069077812 -0.018805924 -0.024590218
## Gpr165   -0.02218366  0.0192150939  0.0041709348 -0.013638851 -0.015836455
## Cbfa2t3  -0.02218013  0.0135865919  0.0076007074  0.010575410 -0.017679377
## Tnfrsf17 -0.02130499  0.0002597742  0.0072914149  0.016976757 -0.006562451
## Il7r     -0.02108119 -0.0049129162  0.0423682248  0.066751928 -0.044075832
## Cryba4   -0.01920738  0.0131492066 -0.0383922518 -0.085652922  0.041492931
## Gbp7     -0.01903686  0.0122138706  0.0314363471 -0.086906707 -0.003969222
## Filip1l  -0.01816183  0.0084317125  0.0288659311  0.009828060 -0.037359629
## Lrba     -0.01709218 -0.0079392318  0.0445312204  0.087127738 -0.062444147
## Ppm1l    -0.01685085 -0.0006411684  0.0120430166  0.021498073 -0.027764018
## Rdx      -0.01563234  0.0322523197 -0.0101063926 -0.028082107 -0.029036583
## Zbtb20   -0.01536199  0.0002726355  0.0438005938  0.048369792 -0.036211015
## Ccr1     -0.01517874  0.0077022669 -0.0120585311 -0.041937733  0.019017410
## Ddit4    -0.01486896 -0.0009223939  0.0121797329  0.023814395 -0.020442643
## Nav2     -0.01409866 -0.0142929110  0.0593419775  0.096965281 -0.061809571
## Arid4a   -0.01399211  0.0061103712  0.0168647372 -0.002183105 -0.013359338
## Igf1r    -0.01345657  0.0071848689 -0.0088389407 -0.006032026 -0.032765500
## Fam102b  -0.01311604 -0.0077763073  0.0377571497  0.088769528 -0.083473741
## Tet1     -0.01295294  0.0037873297  0.0004408433 -0.002861765 -0.013703246
## Khdrbs3  -0.01269506  0.0065244511  0.0181147008  0.048626673 -0.051940118
## Tmem52   -0.01262413 -0.0011493721  0.0012116385 -0.023795250  0.014379272
## Gas6     -0.01261235  0.0120974723 -0.0605701313 -0.088664181 -0.006754880
## Pecam1   -0.01244409  0.0009211247  0.0080679071  0.008475348 -0.029843046
## Upk1b    -0.01235602 -0.0033641098  0.0245727728  0.052303857 -0.031521103
## Tmem176a -0.01177680  0.0057389080  0.0020905264 -0.013429253  0.033399460
##                   PC_6         PC_7          PC_8
## Crybb1    0.0238769004  0.011048220 -0.0028997115
## Ccr5      0.0569038358  0.055592046 -0.0139333431
## Sox4      0.0435266299  0.002963676 -0.0356961419
## Fchsd2    0.0294158288  0.089019658 -0.0232120558
## Hpgd      0.0121485849 -0.105016017 -0.0845414809
## Clec4a3   0.0158861141  0.012877214 -0.0242974813
## Bank1     0.0407457605 -0.001907312 -0.0328508027
## Fcgrt     0.0137658557  0.041666066  0.0123113148
## Lst1     -0.0124565559  0.004216013 -0.0123912356
## Klf7      0.0520305912  0.022738785 -0.0090042175
## Mrc1      0.0533247606  0.121430451  0.0255994223
## Igfbp4    0.0439007800  0.021553984  0.0089468423
## Serpinf1  0.0317907236 -0.009130330 -0.0196456924
## Camk2n1   0.0175620237 -0.027500212 -0.0389277243
## Rgs7bp    0.0392545109  0.017083206 -0.0236065790
## Fry       0.0261683839  0.028172838 -0.0207550320
## Gpr165    0.0475979302  0.008874432  0.0004577951
## Cbfa2t3   0.0264010781 -0.009471263 -0.0377540937
## Tnfrsf17  0.0142571530 -0.023715601 -0.0065563969
## Il7r      0.0023577848 -0.067293856 -0.0658356935
## Cryba4   -0.0142049472  0.055312942 -0.0181824270
## Gbp7      0.0514019814  0.019025147 -0.0137293020
## Filip1l   0.0351766859 -0.010086337 -0.0313252475
## Lrba      0.0202761343 -0.079276963 -0.0676261588
## Ppm1l     0.0192811498 -0.008495790 -0.0228033747
## Rdx       0.0118061661  0.013140204 -0.0256400848
## Zbtb20    0.0048243371  0.022296347 -0.0304752079
## Ccr1      0.0099123317  0.024381823  0.0193727422
## Ddit4     0.0137738355 -0.038357488 -0.0214061887
## Nav2     -0.0108633755 -0.027479937 -0.0837771731
## Arid4a    0.0190071055  0.036486963 -0.0422321973
## Igf1r     0.0076440834  0.047229747 -0.0039480577
## Fam102b   0.0115346031 -0.071793430 -0.1078576193
## Tet1      0.0114894028 -0.003000095 -0.0217079305
## Khdrbs3   0.0148321090 -0.073473649 -0.0534153922
## Tmem52   -0.0002480091  0.011375446  0.0103857805
## Gas6      0.0296092199  0.077575629 -0.0062290112
## Pecam1    0.0091984494 -0.009386596 -0.0231653032
## Upk1b     0.0075888593 -0.053299973 -0.0420004699
## Tmem176a -0.0207360998 -0.003450578 -0.0048027751
decide PCs to use
ElbowPlot(GEX.seur,ndims = 50)

PCs <- 1:20
GEX.seur <- FindNeighbors(GEX.seur, dims = PCs, k.param = 20)
## Computing nearest neighbor graph
## Computing SNN
GEX.seur <- FindClusters(GEX.seur, method = 'igraph' ,resolution = 1.2)
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 23692
## Number of edges: 592986
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.7501
## Number of communities: 16
## Elapsed time: 4 seconds

Run UMAP/tSNE

GEX.seur <- RunTSNE(GEX.seur, dims=PCs, complexity = 100)
GEX.seur <- RunUMAP(GEX.seur, dims=PCs, n.neighbors = 20, seed.use = 287)
## Warning: The default method for RunUMAP has changed from calling Python UMAP via reticulate to the R-native UWOT using the cosine metric
## To use Python UMAP via reticulate, set umap.method to 'umap-learn' and metric to 'correlation'
## This message will be shown once per session
## 15:52:47 UMAP embedding parameters a = 0.9922 b = 1.112
## 15:52:47 Read 23692 rows and found 20 numeric columns
## 15:52:47 Using Annoy for neighbor search, n_neighbors = 20
## 15:52:47 Building Annoy index with metric = cosine, n_trees = 50
## 0%   10   20   30   40   50   60   70   80   90   100%
## [----|----|----|----|----|----|----|----|----|----|
## **************************************************|
## 15:52:50 Writing NN index file to temp file C:\Users\Shaorui\AppData\Local\Temp\RtmpIHLR9L\filedccc2e1a3c91
## 15:52:50 Searching Annoy index using 1 thread, search_k = 2000
## 15:52:56 Annoy recall = 100%
## 15:52:57 Commencing smooth kNN distance calibration using 1 thread
## 15:52:58 Initializing from normalized Laplacian + noise
## 15:52:59 Commencing optimization for 200 epochs, with 683800 positive edges
## 15:53:22 Optimization finished
#saveRDS(GEX.seur,"GEX0811.seur.pure_test.rds")
DimPlot(GEX.seur, reduction = "tsne", label = T) + DimPlot(GEX.seur, reduction = "umap", label = T)

FeaturePlot(GEX.seur, reduction = "umap", features = c("nFeature_RNA","nCount_RNA","percent.mt","percent.rb"))

DimPlot(GEX.seur, reduction = "umap", label = T, group.by = "preAnno") +
  DimPlot(GEX.seur, reduction = "umap", label = F, group.by = "cnt", cols = color.cnt)

DimPlot(GEX.seur, reduction = "umap", label = T,label.size = 2.4, group.by = "preAnno1") +
  DimPlot(GEX.seur, reduction = "umap", label = F, group.by = "preAnno2")

GEX.seur$sort_clusters <- factor(as.character(GEX.seur$seurat_clusters),
                                 levels = c(0,1,9,10,6,8,11,14,
                                            2,3,4,7,12,13,5,15))
VlnPlot(GEX.seur, features = c("nFeature_RNA", "nCount_RNA", "percent.mt", "percent.rb"), 
        ncol = 2,
        pt.size = 0.05, group.by = "sort_clusters")

VlnPlot(GEX.seur, features = c("nFeature_RNA", "nCount_RNA", "percent.mt", "percent.rb"), 
        ncol = 2,
        pt.size = 0, group.by = "sort_clusters")

VlnPlot(GEX.seur, features = c("nFeature_RNA", "nCount_RNA", "percent.mt", "percent.rb"), 
        ncol = 2,
        pt.size = 0.05, group.by = "FB.new", cols = color.FBnew)

VlnPlot(GEX.seur, features = c("nFeature_RNA", "nCount_RNA", "percent.mt", "percent.rb"), 
        ncol = 2,
        pt.size = 0, group.by = "FB.new", cols = color.FBnew)

markers.mig <- c("Ptprc",#"Cd3d","Cd3e","Cd19",
                 "Cd74","Lyz2","Ccl4",
                 "Aif1","P2ry12","C1qa","Spp1",
                 "Top2a","Pcna","Mki67","Mcm6",
                 "Cx3cr1","Il4ra","Il13ra1","Spp1-EGFP",
                 "Fabp5","Hmox1","Ms4a7","Cenpa")
FeaturePlot(GEX.seur, 
            features = markers.mig,
            ncol = 4)

DotPlot(GEX.seur, features = c("Cx3cr1","Spp1","Aif1", # Aif1: Iba1
                               "Fcer1g","Il4ra","Il13ra1"),
        group.by = "sort_clusters")

DotPlot(GEX.seur, features = markers.mig, group.by = "sort_clusters",
                         cols = c("midnightblue","darkorange1")) +
                 theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1))#+ scale_y_discrete(limits=rev)

preAnno

# Anno1
GEX.seur$Anno1 <- as.character(GEX.seur$seurat_clusters)


GEX.seur$Anno1[GEX.seur$Anno1 %in% c(0)] <- "MIG.a1"
GEX.seur$Anno1[GEX.seur$Anno1 %in% c(1)] <- "MIG.a2"
GEX.seur$Anno1[GEX.seur$Anno1 %in% c(9)] <- "MIG.a3"
GEX.seur$Anno1[GEX.seur$Anno1 %in% c(10)] <- "MIG.a4"
GEX.seur$Anno1[GEX.seur$Anno1 %in% c(6)] <- "MIG.a5"

GEX.seur$Anno1[GEX.seur$Anno1 %in% c(8)] <- "MIG.a6"
GEX.seur$Anno1[GEX.seur$Anno1 %in% c(11)] <- "MIG.a7"
GEX.seur$Anno1[GEX.seur$Anno1 %in% c(14)] <- "MIG.a8"

GEX.seur$Anno1[GEX.seur$Anno1 %in% c(2)] <- "MIG.b1"
GEX.seur$Anno1[GEX.seur$Anno1 %in% c(3)] <- "MIG.b2"
GEX.seur$Anno1[GEX.seur$Anno1 %in% c(4)] <- "MIG.b3"
GEX.seur$Anno1[GEX.seur$Anno1 %in% c(7)] <- "MIG.b4"
GEX.seur$Anno1[GEX.seur$Anno1 %in% c(12)] <- "MIG.b5"
GEX.seur$Anno1[GEX.seur$Anno1 %in% c(13)] <- "MIG.b6"
GEX.seur$Anno1[GEX.seur$Anno1 %in% c(5)] <- "MIG.b7"
GEX.seur$Anno1[GEX.seur$Anno1 %in% c(15)] <- "MIG.b8"


GEX.seur$Anno1 <- factor(GEX.seur$Anno1,
                            levels = c(paste0("MIG.a",1:8),
                                       paste0("MIG.b",1:8)))


# Anno2
GEX.seur$Anno2 <- as.character(GEX.seur$Anno1)

GEX.seur$Anno2 <- gsub("1|2|3|4|5|6|7|8","",GEX.seur$Anno2)
GEX.seur$Anno2 <- factor(GEX.seur$Anno2,
                            levels = c("MIG.a","MIG.b"))
pp.new.1a <- DimPlot(GEX.seur, reduction = "umap", label = T, group.by = "seurat_clusters") +
  DimPlot(GEX.seur, reduction = "umap", label = F, group.by = "cnt", cols = color.cnt)
pp.new.1a

pp.new.1b <- DimPlot(GEX.seur, reduction = "umap", label = T,label.size = 3.45, group.by = "Anno1") +
  DimPlot(GEX.seur, reduction = "umap", label = F, group.by = "Anno2", cols = color.cnt[c(3,6)])
pp.new.1b

pp.new.1c <- DimPlot(GEX.seur, reduction = "umap", group.by = "FB.new", split.by = "FB.new", cols = color.FBnew, ncol = 5)
pp.new.1c 

pp.new.1d <- DimPlot(GEX.seur, reduction = "umap", group.by = "cnt", split.by = "cnt", cols = color.cnt, ncol = 3)
pp.new.1d 

markers

# find markers for every cluster compared to all remaining cells, report only the positive ones
Idents(GEX.seur) <- "Anno1"

#GEX.markers.pre <- FindAllMarkers(GEX.seur, only.pos = TRUE, 
#                                  min.pct = 0.05,
#                                  test.use = "MAST",
#                                  logfc.threshold = 0.25)
GEX.markers.pre <- read.table("sc10x_LYN.0811.marker.Anno1.csv", header = TRUE, sep = ",")
GEX.markers.pre %>% group_by(cluster) %>% top_n(n = 8, wt = avg_log2FC)
## # A tibble: 128 x 7
## # Groups:   cluster [16]
##        p_val avg_log2FC pct.1 pct.2 p_val_adj cluster gene         
##        <dbl>      <dbl> <dbl> <dbl>     <dbl> <chr>   <chr>        
##  1 0              0.546 0.994 0.962 0         MIG.a1  Rps19        
##  2 0              0.526 0.997 0.988 0         MIG.a1  Rps11        
##  3 0              0.524 0.999 0.993 0         MIG.a1  Rps27a       
##  4 3.01e-318      0.519 0.994 0.97  5.53e-314 MIG.a1  Rps18        
##  5 1.52e-264      0.558 0.962 0.89  2.79e-260 MIG.a1  Ltc4s        
##  6 5.36e-168      0.591 0.812 0.675 9.85e-164 MIG.a1  Crybb1       
##  7 4.80e-155      0.566 0.657 0.516 8.81e-151 MIG.a1  2410006H16Rik
##  8 2.39e- 95      0.587 0.44  0.29  4.39e- 91 MIG.a1  Mrc1         
##  9 1.18e-273      0.718 0.891 0.706 2.17e-269 MIG.a2  Ccr5         
## 10 6.88e-264      0.720 0.795 0.527 1.26e-259 MIG.a2  Fchsd2       
## # ... with 118 more rows
GEX.markers.anno <- GEX.markers.pre
GEX.markers.anno$cluster <- factor(as.character(GEX.markers.anno$cluster),
                          levels = levels(GEX.seur$Anno1))


markers.pre_t48 <- (GEX.markers.anno %>% group_by(cluster) %>% 
                  filter(pct.1>0.1) %>%
                   top_n(n = 48, wt = avg_log2FC) %>%
                   ungroup() %>%
  arrange(desc(avg_log2FC*pct.1),gene) %>%
                             distinct(gene, .keep_all = TRUE) %>%
                             arrange(cluster,p_val_adj))$gene

markers.pre_t60 <- (GEX.markers.anno %>% group_by(cluster) %>% 
                  filter(pct.1>0.1 & gene %in% grep("Rps|Rpl|Gm|Rik$",GEX.markers.anno$gene,invert = T,value = T)) %>%
                   top_n(n = 60, wt = avg_log2FC) %>%
                   ungroup() %>%
  arrange(desc(avg_log2FC*pct.1),gene) %>%
                             distinct(gene, .keep_all = TRUE) %>%
                             arrange(cluster,p_val_adj))$gene

markers.pre_t120 <- (GEX.markers.anno %>% group_by(cluster) %>% 
                  filter(pct.1>0.1 & gene %in% grep("Rps|Rpl|Gm|Rik$",GEX.markers.anno$gene,invert = T,value = T)) %>%
                   top_n(n = 120, wt = avg_log2FC) %>%
                   ungroup() %>%
  arrange(desc(avg_log2FC*pct.1),gene) %>%
                             distinct(gene, .keep_all = TRUE) %>%
                             arrange(cluster,p_val_adj))$gene
DotPlot(GEX.seur, group.by = "Anno1", features = rev(markers.pre_t60[1:64]))  + coord_flip() + 
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6))

DotPlot(GEX.seur, group.by = "Anno1", features = rev(markers.pre_t60[65:128]))  + coord_flip() + 
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6))

DotPlot(GEX.seur, group.by = "Anno1", features = rev(markers.pre_t60[129:192]))  + coord_flip() + 
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6))

DotPlot(GEX.seur, group.by = "Anno1", features = rev(markers.pre_t60[193:256]))  + coord_flip() + 
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6))

DotPlot(GEX.seur, group.by = "Anno1", features = rev(markers.pre_t60[257:320]))  + coord_flip() + 
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6))

DotPlot(GEX.seur, group.by = "Anno1", features = rev(markers.pre_t60[321:384]))  + coord_flip() + 
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6))

DotPlot(GEX.seur, group.by = "Anno1", features = rev(markers.pre_t60[385:448]))  + coord_flip() + 
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6))

DotPlot(GEX.seur, group.by = "Anno1", features = rev(markers.pre_t60[449:512]))  + coord_flip() + 
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6))

DotPlot(GEX.seur, group.by = "Anno1", features = rev(markers.pre_t60[513:568]))  + coord_flip() + 
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6))

cell composition

usual

cowplot::plot_grid(
pheatmap::pheatmap(table(FB.info=GEX.seur$FB.new,
      clusters=GEX.seur$Anno1),
                   main = "Cell Count",
      gaps_row = c(1,3,5,6,8),
      gaps_col = c(8),
                   cluster_rows = F, cluster_cols = F, display_numbers = T, number_format = "%.0f", legend = F, silent = T)$gtable,

pheatmap::pheatmap(100*rowRatio(table(FB.info=GEX.seur$FB.new,
      clusters=GEX.seur$Anno1)),
                   main = "Cell Ratio",
      gaps_row = c(1,3,5,6,8),
      gaps_col = c(8),
                   cluster_rows = F, cluster_cols = F, display_numbers = T, number_format = "%.1f", legend = F, silent = T)$gtable,
ncol = 1)

glm.nb-anova.LRT

levels(GEX.seur$FB.new)
##  [1] "P05.CTR.1" "P05.MIG.1" "P05.MIG.2" "P05.GFP.1" "P05.GFP.2" "P28.CTR.1"
##  [7] "P28.MIG.1" "P28.MIG.2" "P28.GFP.1" "P28.GFP.2"
GEX.seur$rep <- gsub("P05.CTR.|P05.MIG.|P05.GFP.|P28.CTR.|P28.MIG.|P28.GFP.","rep",as.character(GEX.seur$FB.new))
head(GEX.seur$rep)
## AAACCCAAGAAATCCA-1 AAACCCAAGCGATCGA-1 AAACCCAAGGTAGCCA-1 AAACCCAAGGTTGAGC-1 
##             "rep1"             "rep2"             "rep1"             "rep1" 
## AAACCCAAGTCACTCA-1 AAACCCAAGTGGAATT-1 
##             "rep1"             "rep2"
stat_sort <- GEX.seur@meta.data[,c("cnt","rep","Anno1")]
dim(stat_sort )
## [1] 23692     3
head(stat_sort )
##                        cnt  rep  Anno1
## AAACCCAAGAAATCCA-1 P05.GFP rep1 MIG.a2
## AAACCCAAGCGATCGA-1 P05.MIG rep2 MIG.a2
## AAACCCAAGGTAGCCA-1 P28.GFP rep1 MIG.b7
## AAACCCAAGGTTGAGC-1 P28.GFP rep1 MIG.b3
## AAACCCAAGTCACTCA-1 P28.GFP rep1 MIG.b4
## AAACCCAAGTGGAATT-1 P05.MIG rep2 MIG.a4
stat_sort.s <- stat_sort %>%
  group_by(cnt,rep,Anno1) %>%
  summarise(count=n()) %>%
  mutate(prop= count/sum(count)) %>%
  ungroup()
## `summarise()` has grouped output by 'cnt', 'rep'. You can override using the
## `.groups` argument.
dim(stat_sort.s)
## [1] 151   5
#stat_sort.s
#write.csv(stat_sort.s, "figures0921/2.pure/composition/composition_Anno1.stat.csv")
#write.csv(stat_sort.s %>% reshape2::recast(cnt + rep ~ Anno1, measure.var = 4), 
#          "figures0921/2.pure/composition/composition_Anno1.stat_df.count.csv")
#write.csv(stat_sort.s %>% reshape2::recast(cnt + rep ~ Anno1, measure.var = 5), 
#          "figures0921/2.pure/composition/composition_Anno1.stat_df.ratio.csv")

plot

pcom.sort <- stat_sort.s %>% #filter(stage %in% c("P00","P07","P14","P21","P28")) %>%
  ggplot(aes(x = Anno1, y = 100*prop, color=cnt)) +
  geom_bar(stat="summary", fun="mean", position = position_dodge(0.75), width = 0.58, fill="white") +
  theme_classic(base_size = 15) + 
  scale_color_manual(values = color.cnt, name="") +
  labs(title="Cell Composition", y="Proportion") + 
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6)) +
  
  geom_point(aes(x=Anno1, y=100*prop, color=cnt),
             position = position_dodge(0.75),
             shape=16,alpha=0.75,size=1.75,
             stroke=0.15, show.legend = F)
pcom.sort

significance

some condition no replicates, not run

test some markers

marker.fig0921 <- c("Ptprc","P2ry12","Spp1","Spp1-EGFP",
                    "Cx3cr1","Csf1r","Aif1","C1qa",
                    "Ccl4","Lpl","Fabp5","Csf1",
                    "Ctsf","Ccr5","Arsb","Itgb1",
                    "Cd81","Lag3","Cd34","H2-K1",
                    "Top2a","Pcna","Mki67","Mcm6"
                    #"Cd74","H2-Aa","Ccr2","Cd209a",
                    #"Mrc1","Pf4","Dab2","Lyve1",
                    #"Tubb3","Actl6b","Syt1","Sst"
                    )

pp.new.2a <- FeaturePlot(GEX.seur, features = marker.fig0921, ncol = 4,raster = T, pt.size = 2.75)
pp.new.2a

pp.new.2b <- DotPlot(GEX.seur, features = marker.fig0921, group.by = "Anno1",
                         cols = c("midnightblue","darkorange1")) +
                 theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1))#+ scale_y_discrete(limits=rev)
pp.new.2b

pp.new.2c <- VlnPlot(GEX.seur, features = marker.fig0921, ncol = 4, pt.size = 0, raster = T) #& geom_jitter(alpha=0.15, shape=16, width = 0.2, size = 0.02)
pp.new.2c

pp.new.2c1 <- VlnPlot(GEX.seur, features = marker.fig0921, ncol = 4, pt.size = 0, raster = F) & geom_jitter(alpha=0.15, shape=16, width = 0.2, size = 0.02)
pp.new.2c1

signature score

#### 10x data, calculate signature score       
#
## The code below is from Adam Hamber
## 2D scoring by Itay
get_controls <- function(counts, gene.list, verbose=F, control.genes.per.gene=10)
{
    # Itay: "Such scores are inevitably correlated with cell complexity so to avoid 
    # that I subtract a "control" score which is generated by averaging over a control 
    # gene set. Control gene sets are chosen to contain 100 times more genes than the 
    # real gene set (analogous to averaging over 100 control sets of similar size) and 
    # to have the same distribution of population/bulk - based expression levels as the 
    # real gene set, such that they are expected to have the same number of "zeros" and 
    # to eliminate the correlation with complexity."
    # ---------------------------------------------------------------------------------
    # Going to find control points by finding the closest genes in terms of expression level and % of the time we observe it
    if(verbose){cat(sprintf("Finding %s background genes based on similarity to given gene set [%s genes] \n", 
                            control.genes.per.gene*length(gene.list), length(gene.list)))}
    cat("Summarizing data \n")
    summary = data.frame(gene=row.names(counts), mean.expr = Matrix::rowMeans(counts), fract.zero = Matrix::rowMeans(counts==0), stringsAsFactors = F)
    #summary = data.frame(gene=row.names(counts), mean.expr = apply(counts,1,mean), fract.zero = apply(counts==0,1,mean), stringsAsFactors = F)
    summary$mean.expr.s = scale(summary$mean.expr)
    summary$fract.zero.s = scale(summary$fract.zero)
    actual.genes = summary[summary$gene %in% gene.list,]
    background.genes = summary[!summary$gene %in% gene.list,]
    
    #find the 10 closest genes to each cell cycle marker gene and add them to the lists of control genes
    get_closest_genes <- function(i)
    {
        background.genes$dist = sqrt((background.genes$mean.expr.s - actual.genes$mean.expr.s[i])^2 + 
                                         (background.genes$fract.zero.s - actual.genes$fract.zero.s[i])^2)
        ordered = background.genes$gene[order(background.genes$dist)]
        ordered = ordered[!ordered %in% controls] # don't take genes that already appear in the list 
        closest = head(ordered, n=control.genes.per.gene)
        return(closest)
    }
    controls = c();
    
    for (i in 1:length(gene.list)){
        #info(sprintf("Finding %s control genes for %s", control.genes.per.gene, gene.list[i]))
        closest = get_closest_genes(i)
        #info(sprintf("Found %s: ", length(closest)))
        controls = unique(c(controls, closest))
    }
    
    if(verbose){cat(sprintf("Control gene selection complete. %s genes found. \n", length(controls)))}
    #print(controls)
    return(controls)
}

## Define calculate function
calculate_signature_score <- function(count_matrix, gene_list){
    control_gene <- get_controls(counts = count_matrix,
                                 gene.list = gene_list)
    signature_score <- colMeans(count_matrix[gene_list, ], na.rm = TRUE) - 
        colMeans(count_matrix[control_gene, ], na.rm = TRUE)
    return(signature_score)
}

add_geneset_score <- function(obj, geneset, setname){
  score <- calculate_signature_score(as.data.frame(obj@assays[['RNA']]@data),
                                     geneset)
  obj <- AddMetaData(obj,
                     score,
                     setname)
  return(obj)
}
## proc_DEG()
#       to process edgeR result for DEGs' comparison
#              mat_cut is a matrix after advanced filtering, 'like TPM > n in at least one condition'
# 

proc_DEG <- function(deg, p.cut=0.05, FC.cut = 2, padj=TRUE, abs=TRUE, mat_cut=NULL, gene_cut=NULL){
    rownames(deg) <- deg$gene
    
    if(padj==TRUE){
        deg <- deg %>% filter(padj < p.cut)
    }else{
        deg <- deg %>% filter(pvalue < p.cut)
    }
    
    if(abs==TRUE){
        deg <- deg %>% filter(abs(FC) > FC.cut)
    }else if(FC.cut >0){
        deg <- deg %>% filter(FC > FC.cut)
    }else{
        deg <- deg %>% filter(FC < FC.cut)
    }
    
    if(!is.null(mat_cut)){
        deg <- deg[rownames(deg) %in% rownames(mat_cut),]
    }
    if(!is.null(gene_cut)){
        deg <- deg[rownames(deg) %in% gene_cut,]
    }
    return(deg)
}

SS2 DEGs

DEGs.SS2_P06 <- read.csv("../bulkDEGs/DEG_table.P06_SIMPLE.Spp1_PvsN.csv",row.names = 1)
DEGs.SS2_P06$Gene <- rownames(DEGs.SS2_P06)
#DEGs.SS2_P06
DEGs.SS2_P28 <- read.csv("../bulkDEGs/DEG_table.P28_SIMPLE.Spp1_PvsN.csv",row.names = 1)
DEGs.SS2_P28$Gene <- rownames(DEGs.SS2_P28)
#DEGs.SS2_P28
DEGs.SS2_P06["Spp1",]
##      log2FoldChange   logCPM       LR   pvalue     padj       FC Gene
## Spp1        6.05263 7.836443 232.6566 1.57e-52 7.95e-49 66.37783 Spp1
DEGs.SS2_P28["Spp1",]
##      log2FoldChange   logCPM       LR   pvalue     padj       FC Gene
## Spp1       9.978018 10.00778 146.9156 8.19e-34 6.34e-32 1008.516 Spp1
DEGs.list <-  list(P06.Spp1_up = proc_DEG(DEGs.SS2_P06, p.cut = 0.05, FC.cut = 1.5, padj = T, abs = FALSE)$Gene,
                   P28.Spp1_up = proc_DEG(DEGs.SS2_P28, p.cut = 0.05, FC.cut = 1.5, padj = T, abs = FALSE)$Gene)
lapply(DEGs.list, length)
## $P06.Spp1_up
## [1] 177
## 
## $P28.Spp1_up
## [1] 1259
pp.DEG <- venn::venn(DEGs.list,
           zcolor = 'style',
           ggplot = T) + 
           labs( title=" cutoff: padj<0.05, FC>1.5")+ 
           theme(plot.title = element_text(size=15))
pp.DEG 

DEG.comp <- list(Spp1_up.P06only = setdiff(DEGs.list$P06.Spp1_up, DEGs.list$P28.Spp1_up),
                 Spp1_up.both = intersect(DEGs.list$P06.Spp1_up, DEGs.list$P28.Spp1_up),
                 Spp1_up.P28only = setdiff(DEGs.list$P28.Spp1_up, DEGs.list$P06.Spp1_up))
lapply(DEG.comp, length)
## $Spp1_up.P06only
## [1] 60
## 
## $Spp1_up.both
## [1] 117
## 
## $Spp1_up.P28only
## [1] 1142
names(DEG.comp)
## [1] "Spp1_up.P06only" "Spp1_up.both"    "Spp1_up.P28only"
for(NN in names(DEG.comp)){
  
  write.table(DEG.comp[[NN]],paste0("figures0921/2.pure/signature/overlap/DEGs.",NN,".num",length(DEG.comp[[NN]]),".txt"), col.names = F, row.names = F ,quote = F)
  
}
GEX.seur <- add_geneset_score(GEX.seur, DEG.comp$Spp1_up.P06only, "Spp1_up.P06only")
## Summarizing data
GEX.seur <- add_geneset_score(GEX.seur, DEG.comp$Spp1_up.both, "Spp1_up.both")
## Summarizing data
GEX.seur <- add_geneset_score(GEX.seur, DEG.comp$Spp1_up.P28only, "Spp1_up.P28only")
## Summarizing data

plot

samples

VlnPlot(GEX.seur, features = c("Spp1_up.P06only","Spp1_up.both","Spp1_up.P28only"), group.by = "FB.new", cols = color.FBnew, ncol = 3)

VlnPlot(GEX.seur, features = c("Spp1_up.P06only","Spp1_up.both","Spp1_up.P28only"), group.by = "cnt", cols = color.cnt, ncol = 3, pt.size = 0) & 
    geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
  stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) &
  ggpubr::stat_compare_means(aes(lable = ..p.signif..), 
                               method = "wilcox.test",
                               comparisons = list(c("P05.MIG","P05.GFP"),
                                                  c("P28.MIG","P28.GFP")),
                               label.y = c(0.3,0.3),
                               size=3.5
                               )

pcomp.1a <- VlnPlot(GEX.seur, features = c("Spp1_up.P06only","Spp1_up.both","Spp1_up.P28only")[1], group.by = "cnt", cols = color.cnt, ncol = 1, pt.size = 0)  + NoLegend() + 
  labs(x="",y="Sigature Score of DEGs") & 
    geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
  stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) &
  ggpubr::stat_compare_means(aes(lable = ..p.signif..), 
                               method = "wilcox.test",
                               comparisons = list(c("P05.MIG","P05.GFP"),
                                                  c("P05.CTR","P05.GFP"),
                                                  c("P28.MIG","P28.GFP"),
                                                  c("P28.CTR","P28.GFP")),
                               label.y = c(0.25,0.3,0.25,0.3),
                               size=2
                               )
pcomp.1a

pcomp.1b <- VlnPlot(GEX.seur, features = c("Spp1_up.P06only","Spp1_up.both","Spp1_up.P28only")[2], group.by = "cnt", cols = color.cnt, ncol = 1, pt.size = 0) + NoLegend() + 
  labs(x="",y="Sigature Score of DEGs") & 
    geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
  stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) &
  ggpubr::stat_compare_means(aes(lable = ..p.signif..), 
                               method = "wilcox.test",
                               comparisons = list(c("P05.MIG","P05.GFP"),
                                                  c("P05.CTR","P05.GFP"),
                                                  c("P28.MIG","P28.GFP"),
                                                  c("P28.CTR","P28.GFP")),
                               label.y = c(0.35,0.5,0.45,0.55),
                               size=2
                               )
pcomp.1b

pcomp.1c <- VlnPlot(GEX.seur, features = c("Spp1_up.P06only","Spp1_up.both","Spp1_up.P28only")[3], group.by = "cnt", cols = color.cnt, ncol = 1, pt.size = 0) + NoLegend() + 
  labs(x="",y="Sigature Score of DEGs") & 
    geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
  stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) &
  ggpubr::stat_compare_means(aes(lable = ..p.signif..), 
                               method = "wilcox.test",
                               comparisons = list(c("P05.MIG","P05.GFP"),
                                                  c("P05.CTR","P05.GFP"),
                                                  c("P28.MIG","P28.GFP"),
                                                  c("P28.CTR","P28.GFP")),
                               label.y = c(0.5,0.55,0.55,0.6),
                               size=2
                               )
pcomp.1c

#
ggsave("figures0921/2.pure/signature/sig_score.DEGs.samples.Spp1_up.P06only.pdf",
       width = 3.5,height = 4,
       plot = pcomp.1a)
ggsave("figures0921/2.pure/signature/sig_score.DEGs.samples.Spp1_up.both.pdf",
       width = 3.5,height = 4,
       plot = pcomp.1b)
ggsave("figures0921/2.pure/signature/sig_score.DEGs.samples.Spp1_up.P28only.pdf",
       width = 3.5,height = 4,
       plot = pcomp.1c)

Anno1

VlnPlot(GEX.seur, features = c("Spp1_up.P06only","Spp1_up.both","Spp1_up.P28only"), group.by = "Anno1", ncol = 3)

VlnPlot(GEX.seur, features = c("Spp1_up.P06only","Spp1_up.both","Spp1_up.P28only"), group.by = "Anno1",  ncol = 3, pt.size = 0) & 
    geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
  stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) &
  ggpubr::stat_compare_means(aes(lable = ..p.signif..), 
                               method = "wilcox.test",
                               comparisons = list(c("P05.MIG","P05.GFP"),
                                                  c("P28.MIG","P28.GFP")),
                               label.y = c(0.3,0.3),
                               size=3.5
                               )
## Warning: Computation failed in `stat_signif()`:
## 需要TRUE/FALSE值的地方不可以用缺少值

## Warning: Computation failed in `stat_signif()`:
## 需要TRUE/FALSE值的地方不可以用缺少值

## Warning: Computation failed in `stat_signif()`:
## 需要TRUE/FALSE值的地方不可以用缺少值

pcomp.2a <- VlnPlot(GEX.seur, features = c("Spp1_up.P06only","Spp1_up.both","Spp1_up.P28only")[1], group.by = "Anno1", ncol = 1, pt.size = 0)  + NoLegend() + 
  labs(x="",y="Sigature Score of DEGs") & 
    geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
  stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) &
  ggpubr::stat_compare_means(aes(lable = ..p.signif..), 
                               method = "wilcox.test",
                               comparisons = list(c("MIG.a3","MIG.a5"),
                                                  c("MIG.a4","MIG.a5"),
                                                  c("MIG.a6","MIG.a5"),
                                                  c("MIG.b5","MIG.b6"),
                                                  c("MIG.b4","MIG.b6"),
                                                  c("MIG.b6","MIG.b7"),
                                                  c("MIG.b6","MIG.b8")),
                               label.y = c(0.25,0.2,0.17,0.25,0.32,0.23,0.28),
                               size=2
                               )
pcomp.2a

pcomp.2b <- VlnPlot(GEX.seur, features = c("Spp1_up.P06only","Spp1_up.both","Spp1_up.P28only")[2], group.by = "Anno1", ncol = 1, pt.size = 0) + NoLegend() + 
  labs(x="",y="Sigature Score of DEGs") & 
    geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
  stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) &
  ggpubr::stat_compare_means(aes(lable = ..p.signif..), 
                               method = "wilcox.test",
                               comparisons = list(c("MIG.a3","MIG.a5"),
                                                  c("MIG.a4","MIG.a5"),
                                                  c("MIG.a6","MIG.a5"),
                                                  c("MIG.b5","MIG.b6"),
                                                  c("MIG.b4","MIG.b6"),
                                                  c("MIG.b6","MIG.b7"),
                                                  c("MIG.b6","MIG.b8")),
                               label.y = c(0.4,0.35,0.3,0.48,0.55,0.52,0.6),
                               size=2
                               )
pcomp.2b

pcomp.2c <- VlnPlot(GEX.seur, features = c("Spp1_up.P06only","Spp1_up.both","Spp1_up.P28only")[3], group.by = "Anno1", ncol = 1, pt.size = 0) + NoLegend() + 
  labs(x="",y="Sigature Score of DEGs") & 
    geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
  stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) &
  ggpubr::stat_compare_means(aes(lable = ..p.signif..), 
                               method = "wilcox.test",
                               comparisons = list(c("MIG.a3","MIG.a5"),
                                                  c("MIG.a4","MIG.a5"),
                                                  c("MIG.a6","MIG.a5"),
                                                  c("MIG.b5","MIG.b6"),
                                                  c("MIG.b4","MIG.b6"),
                                                  c("MIG.b6","MIG.b7"),
                                                  c("MIG.b6","MIG.b8")),
                               label.y = c(0.55,0.5,0.47,0.575,0.6,0.55,0.59),
                               size=2
                               )
pcomp.2c

#
ggsave("figures0921/2.pure/signature/sig_score.DEGs.Anno1.Spp1_up.P06only.pdf",
       width = 7.5,height = 4,
       plot = pcomp.2a)
ggsave("figures0921/2.pure/signature/sig_score.DEGs.Anno1.Spp1_up.both.pdf",
       width = 7.5,height = 4,
       plot = pcomp.2b)
ggsave("figures0921/2.pure/signature/sig_score.DEGs.Anno1.Spp1_up.P28only.pdf",
       width = 7.5,height = 4,
       plot = pcomp.2c)
#saveRDS(GEX.seur,"./GEX0811.seur.pure_Anno1.rds")

forVelocity

head(GEX.seur,4)
##                     orig.ident nCount_RNA nFeature_RNA percent.mt percent.rb
## AAACCCAAGAAATCCA-1 plus_SIMPLE       2998         1458   2.701801   15.01001
## AAACCCAAGCGATCGA-1 plus_SIMPLE       5474         2107   5.608330   17.97589
## AAACCCAAGGTAGCCA-1 plus_SIMPLE       7120         2540   3.609551   12.48596
## AAACCCAAGGTTGAGC-1 plus_SIMPLE       4650         1705   4.430108   19.03226
##                      FB.info     S.Score   G2M.Score Phase RNA_snn_res.1.5
## AAACCCAAGAAATCCA-1 P05.GFP.1 -0.06682298 -0.09221214    G1               2
## AAACCCAAGCGATCGA-1 P05.MIG.2 -0.09983915 -0.03784219    G1               1
## AAACCCAAGGTAGCCA-1 P28.GFP.1 -0.07149868 -0.02169489    G1               9
## AAACCCAAGGTTGAGC-1 P28.GFP.1  0.03069876 -0.11000513     S               5
##                    seurat_clusters RNA_snn_res.1.2     cnt sort_clusters
## AAACCCAAGAAATCCA-1               1               1 P05.GFP             1
## AAACCCAAGCGATCGA-1               1               1 P05.MIG             1
## AAACCCAAGGTAGCCA-1               5               5 P28.GFP             5
## AAACCCAAGGTTGAGC-1               4               4 P28.GFP             4
##                    pANN_0.25_0.005_1267 DoubletFinder0.05 pANN_0.25_0.005_2534
## AAACCCAAGAAATCCA-1            0.0000000           Singlet            0.0000000
## AAACCCAAGCGATCGA-1            0.4970414           Singlet            0.4911243
## AAACCCAAGGTAGCCA-1            0.4437870           Singlet            0.4260355
## AAACCCAAGGTTGAGC-1            0.1420118           Singlet            0.1183432
##                    DoubletFinder0.1 preAnno1 preAnno2 preAnno    FB.new  Anno1
## AAACCCAAGAAATCCA-1          Singlet   MIG.a2    MIG.a     MIG P05.GFP.1 MIG.a2
## AAACCCAAGCGATCGA-1          Doublet   MIG.a2    MIG.a     MIG P05.MIG.2 MIG.a2
## AAACCCAAGGTAGCCA-1          Doublet   MIG.b7    MIG.b     MIG P28.GFP.1 MIG.b7
## AAACCCAAGGTTGAGC-1          Singlet   MIG.b7    MIG.b     MIG P28.GFP.1 MIG.b3
##                    Anno2  rep Spp1_up.P06only Spp1_up.both Spp1_up.P28only
## AAACCCAAGAAATCCA-1 MIG.a rep1      0.07799678   0.01572847       0.2361100
## AAACCCAAGCGATCGA-1 MIG.a rep2     -0.04681886  -0.05326058       0.2947471
## AAACCCAAGGTAGCCA-1 MIG.b rep1     -0.03821920   0.29210607       0.4436883
## AAACCCAAGGTTGAGC-1 MIG.b rep1     -0.00965059   0.10195132       0.4066607

output

mainly follow the official scVelo tutorial
https://smorabit.github.io/tutorials/8_velocyto/

save metadata table

## re-process cell barcodes so that next would be easier in jupyter notebook
#GEX.seur$tissue <- sapply(GEX.seur$orig.ident, function(x){base::strsplit(x,split=".",fixed = T)[[1]][2]})

#GEX.seur$barcode <- paste0(sapply(rownames(GEX.seur@meta.data), function(x){base::strsplit(x,split="-",fixed = T)[[1]][1]}),
#                           ".",
#                           GEX.seur$tissue)

GEX.seur$barcode <- rownames(GEX.seur@meta.data)

GEX.seur$UMAP_1 <- GEX.seur@reductions$umap@cell.embeddings[,1]
GEX.seur$UMAP_2 <- GEX.seur@reductions$umap@cell.embeddings[,2]

head(GEX.seur@meta.data)
##                     orig.ident nCount_RNA nFeature_RNA percent.mt percent.rb
## AAACCCAAGAAATCCA-1 plus_SIMPLE       2998         1458   2.701801   15.01001
## AAACCCAAGCGATCGA-1 plus_SIMPLE       5474         2107   5.608330   17.97589
## AAACCCAAGGTAGCCA-1 plus_SIMPLE       7120         2540   3.609551   12.48596
## AAACCCAAGGTTGAGC-1 plus_SIMPLE       4650         1705   4.430108   19.03226
## AAACCCAAGTCACTCA-1 plus_SIMPLE       2150          894   2.790698   10.04651
## AAACCCAAGTGGAATT-1 plus_SIMPLE       3738         1568   7.276619   20.43874
##                      FB.info     S.Score   G2M.Score Phase RNA_snn_res.1.5
## AAACCCAAGAAATCCA-1 P05.GFP.1 -0.06682298 -0.09221214    G1               2
## AAACCCAAGCGATCGA-1 P05.MIG.2 -0.09983915 -0.03784219    G1               1
## AAACCCAAGGTAGCCA-1 P28.GFP.1 -0.07149868 -0.02169489    G1               9
## AAACCCAAGGTTGAGC-1 P28.GFP.1  0.03069876 -0.11000513     S               5
## AAACCCAAGTCACTCA-1 P28.GFP.1 -0.01015567 -0.02782644    G1              11
## AAACCCAAGTGGAATT-1 P05.MIG.2 -0.02148368 -0.04812337    G1               6
##                    seurat_clusters RNA_snn_res.1.2     cnt sort_clusters
## AAACCCAAGAAATCCA-1               1               1 P05.GFP             1
## AAACCCAAGCGATCGA-1               1               1 P05.MIG             1
## AAACCCAAGGTAGCCA-1               5               5 P28.GFP             5
## AAACCCAAGGTTGAGC-1               4               4 P28.GFP             4
## AAACCCAAGTCACTCA-1               7               7 P28.GFP             7
## AAACCCAAGTGGAATT-1              10              10 P05.MIG            10
##                    pANN_0.25_0.005_1267 DoubletFinder0.05 pANN_0.25_0.005_2534
## AAACCCAAGAAATCCA-1           0.00000000           Singlet           0.00000000
## AAACCCAAGCGATCGA-1           0.49704142           Singlet           0.49112426
## AAACCCAAGGTAGCCA-1           0.44378698           Singlet           0.42603550
## AAACCCAAGGTTGAGC-1           0.14201183           Singlet           0.11834320
## AAACCCAAGTCACTCA-1           0.00000000           Singlet           0.00000000
## AAACCCAAGTGGAATT-1           0.01183432           Singlet           0.01775148
##                    DoubletFinder0.1 preAnno1 preAnno2 preAnno    FB.new  Anno1
## AAACCCAAGAAATCCA-1          Singlet   MIG.a2    MIG.a     MIG P05.GFP.1 MIG.a2
## AAACCCAAGCGATCGA-1          Doublet   MIG.a2    MIG.a     MIG P05.MIG.2 MIG.a2
## AAACCCAAGGTAGCCA-1          Doublet   MIG.b7    MIG.b     MIG P28.GFP.1 MIG.b7
## AAACCCAAGGTTGAGC-1          Singlet   MIG.b7    MIG.b     MIG P28.GFP.1 MIG.b3
## AAACCCAAGTCACTCA-1          Singlet   MIG.b6    MIG.b     MIG P28.GFP.1 MIG.b4
## AAACCCAAGTGGAATT-1          Singlet   MIG.a3    MIG.a     MIG P05.MIG.2 MIG.a4
##                    Anno2  rep Spp1_up.P06only Spp1_up.both Spp1_up.P28only
## AAACCCAAGAAATCCA-1 MIG.a rep1      0.07799678   0.01572847       0.2361100
## AAACCCAAGCGATCGA-1 MIG.a rep2     -0.04681886  -0.05326058       0.2947471
## AAACCCAAGGTAGCCA-1 MIG.b rep1     -0.03821920   0.29210607       0.4436883
## AAACCCAAGGTTGAGC-1 MIG.b rep1     -0.00965059   0.10195132       0.4066607
## AAACCCAAGTCACTCA-1 MIG.b rep1      0.03137777   0.29949531       0.3521602
## AAACCCAAGTGGAATT-1 MIG.a rep2     -0.02037658  -0.01007519       0.3301183
##                               barcode     UMAP_1    UMAP_2
## AAACCCAAGAAATCCA-1 AAACCCAAGAAATCCA-1  1.1502470  4.239640
## AAACCCAAGCGATCGA-1 AAACCCAAGCGATCGA-1 -0.9110060  3.796429
## AAACCCAAGGTAGCCA-1 AAACCCAAGGTAGCCA-1 -0.1814226 -7.373271
## AAACCCAAGGTTGAGC-1 AAACCCAAGGTTGAGC-1  0.4916820 -4.882718
## AAACCCAAGTCACTCA-1 AAACCCAAGTCACTCA-1 -0.8550854 -3.813902
## AAACCCAAGTGGAATT-1 AAACCCAAGTGGAATT-1 -1.9329503  3.032112
sum(duplicated(GEX.seur$barcode))
## [1] 0

write counts matrix

write pca

write filtered genes

figures1002

GEX.seur <- readRDS("./GEX0811.seur.pure_Anno1.rds")
GEX.seur
## An object of class Seurat 
## 18371 features across 23692 samples within 1 assay 
## Active assay: RNA (18371 features, 1389 variable features)
##  3 dimensional reductions calculated: pca, tsne, umap
head(GEX.seur@meta.data)
##                     orig.ident nCount_RNA nFeature_RNA percent.mt percent.rb
## AAACCCAAGAAATCCA-1 plus_SIMPLE       2998         1458   2.701801   15.01001
## AAACCCAAGCGATCGA-1 plus_SIMPLE       5474         2107   5.608330   17.97589
## AAACCCAAGGTAGCCA-1 plus_SIMPLE       7120         2540   3.609551   12.48596
## AAACCCAAGGTTGAGC-1 plus_SIMPLE       4650         1705   4.430108   19.03226
## AAACCCAAGTCACTCA-1 plus_SIMPLE       2150          894   2.790698   10.04651
## AAACCCAAGTGGAATT-1 plus_SIMPLE       3738         1568   7.276619   20.43874
##                      FB.info     S.Score   G2M.Score Phase RNA_snn_res.1.5
## AAACCCAAGAAATCCA-1 P05.GFP.1 -0.06682298 -0.09221214    G1               2
## AAACCCAAGCGATCGA-1 P05.MIG.2 -0.09983915 -0.03784219    G1               1
## AAACCCAAGGTAGCCA-1 P28.GFP.1 -0.07149868 -0.02169489    G1               9
## AAACCCAAGGTTGAGC-1 P28.GFP.1  0.03069876 -0.11000513     S               5
## AAACCCAAGTCACTCA-1 P28.GFP.1 -0.01015567 -0.02782644    G1              11
## AAACCCAAGTGGAATT-1 P05.MIG.2 -0.02148368 -0.04812337    G1               6
##                    seurat_clusters RNA_snn_res.1.2     cnt sort_clusters
## AAACCCAAGAAATCCA-1               1               1 P05.GFP             1
## AAACCCAAGCGATCGA-1               1               1 P05.MIG             1
## AAACCCAAGGTAGCCA-1               5               5 P28.GFP             5
## AAACCCAAGGTTGAGC-1               4               4 P28.GFP             4
## AAACCCAAGTCACTCA-1               7               7 P28.GFP             7
## AAACCCAAGTGGAATT-1              10              10 P05.MIG            10
##                    pANN_0.25_0.005_1267 DoubletFinder0.05 pANN_0.25_0.005_2534
## AAACCCAAGAAATCCA-1           0.00000000           Singlet           0.00000000
## AAACCCAAGCGATCGA-1           0.49704142           Singlet           0.49112426
## AAACCCAAGGTAGCCA-1           0.44378698           Singlet           0.42603550
## AAACCCAAGGTTGAGC-1           0.14201183           Singlet           0.11834320
## AAACCCAAGTCACTCA-1           0.00000000           Singlet           0.00000000
## AAACCCAAGTGGAATT-1           0.01183432           Singlet           0.01775148
##                    DoubletFinder0.1 preAnno1 preAnno2 preAnno    FB.new  Anno1
## AAACCCAAGAAATCCA-1          Singlet   MIG.a2    MIG.a     MIG P05.GFP.1 MIG.a2
## AAACCCAAGCGATCGA-1          Doublet   MIG.a2    MIG.a     MIG P05.MIG.2 MIG.a2
## AAACCCAAGGTAGCCA-1          Doublet   MIG.b7    MIG.b     MIG P28.GFP.1 MIG.b7
## AAACCCAAGGTTGAGC-1          Singlet   MIG.b7    MIG.b     MIG P28.GFP.1 MIG.b3
## AAACCCAAGTCACTCA-1          Singlet   MIG.b6    MIG.b     MIG P28.GFP.1 MIG.b4
## AAACCCAAGTGGAATT-1          Singlet   MIG.a3    MIG.a     MIG P05.MIG.2 MIG.a4
##                    Anno2  rep Spp1_up.P06only Spp1_up.both Spp1_up.P28only
## AAACCCAAGAAATCCA-1 MIG.a rep1      0.07799678   0.01572847       0.2361100
## AAACCCAAGCGATCGA-1 MIG.a rep2     -0.04681886  -0.05326058       0.2947471
## AAACCCAAGGTAGCCA-1 MIG.b rep1     -0.03821920   0.29210607       0.4436883
## AAACCCAAGGTTGAGC-1 MIG.b rep1     -0.00965059   0.10195132       0.4066607
## AAACCCAAGTCACTCA-1 MIG.b rep1      0.03137777   0.29949531       0.3521602
## AAACCCAAGTGGAATT-1 MIG.a rep2     -0.02037658  -0.01007519       0.3301183
##                    DAM.sig_top50 DAM.sig_top100 DAM.sig_top250 DAM.sig_top500
## AAACCCAAGAAATCCA-1   -0.15072706    -0.18592966    -0.13817389     0.02798874
## AAACCCAAGCGATCGA-1   -0.23592021    -0.16868345    -0.06001030     0.12662573
## AAACCCAAGGTAGCCA-1    0.79968333     0.61875031     0.48767439     0.50833408
## AAACCCAAGGTTGAGC-1    0.25948931     0.20225591     0.17024206     0.31339640
## AAACCCAAGTCACTCA-1    0.90775071     0.64612166     0.38459858     0.38874980
## AAACCCAAGTGGAATT-1   -0.02039969    -0.04251335    -0.02195155     0.17200874
##                    Cycling_score
## AAACCCAAGAAATCCA-1  -0.072838959
## AAACCCAAGCGATCGA-1  -0.045974881
## AAACCCAAGGTAGCCA-1  -0.006387924
## AAACCCAAGGTTGAGC-1  -0.010882599
## AAACCCAAGTCACTCA-1  -0.036946642
## AAACCCAAGTGGAATT-1  -0.007985408

ppnew.1

Signature plot as feature plot

ppnew.1a <- FeaturePlot(GEX.seur,features = c("Spp1_up.P06only","Spp1_up.both","Spp1_up.P28only"), ncol = 3,
                        raster = T, pt.size = 3.5,
            min.cutoff = -0.18, max.cutoff = 0.8, keep.scale = "all")
ppnew.1a

ppnew.1b <- FeaturePlot(GEX.seur,features = c("Spp1_up.P06only","Spp1_up.both","Spp1_up.P28only"), ncol = 3,
                        cols = c("lightgrey","red"),raster = T, pt.size = 3.5,
            min.cutoff = -0.18, max.cutoff = 0.8, keep.scale = "all")
ppnew.1b

mapal <- colorRampPalette(RColorBrewer::brewer.pal(4,"Spectral"))(120)
ppnew.1c <- FeaturePlot(GEX.seur,features = c("Spp1_up.P06only","Spp1_up.both","Spp1_up.P28only"), ncol = 3,
                        #cols = c("lightgrey","red"), 
                        #cols = rev(mapal),
                        raster = T, pt.size = 3.5,
            min.cutoff = -0.18, max.cutoff = 0.8, keep.scale = "all") & scale_color_gradientn(colors = rev(mapal), limits= c(-0.16,0.65), breaks = c(0.0,0.2,0.4,0.6))
ppnew.1c

ppnew.1d <- FeaturePlot(GEX.seur,features = c("Spp1_up.P06only","Spp1_up.both","Spp1_up.P28only"), ncol = 3,
                        #cols = c("lightgrey","red"), 
                        #cols = rev(mapal),
                        raster = T, pt.size = 3.5,
            min.cutoff = -0.18, max.cutoff = 0.8, keep.scale = "all") & scale_color_gradientn(colors = rev(mapal))
ppnew.1d

mod

using own scale, highlight the hot-regions

ppnew.1d1 <- FeaturePlot(GEX.seur,features = c("Spp1_up.P06only"), ncol = 1,
                        #cols = c("lightgrey","red"), 
                        #cols = rev(mapal),
                        raster = T, pt.size = 3.5,
            min.cutoff = -0.08, max.cutoff = 0.38) & scale_color_gradientn(colors = rev(mapal), #limits= c(-0.08,0.38), 
                                                                           breaks=c(0.0,0.1,0.2,0.3))
ppnew.1d1

ppnew.1d2 <- FeaturePlot(GEX.seur,features = c("Spp1_up.both"), ncol = 1,
                        #cols = c("lightgrey","red"), 
                        #cols = rev(mapal),
                        raster = T, pt.size = 3.5,
            min.cutoff = -0.16, max.cutoff = 0.46) + scale_color_gradientn(colors = rev(mapal), #limits= c(-0.16,0.46), 
                                                                                              breaks=c(0.0,0.2,0.4))
ppnew.1d2

ppnew.1d3a <- FeaturePlot(GEX.seur,features = c("Spp1_up.P28only"), ncol = 1,
                        #cols = c("lightgrey","red"), 
                        #cols = rev(mapal),
                        raster = T, pt.size = 3.5,
            min.cutoff = -0.1, max.cutoff = 0.54) & scale_color_gradientn(colors = rev(mapal), #limits= c(-0.1,0.64), 
                                                                          breaks=c(0.1,0.2,0.3,0.4,0.5))
ppnew.1d3a

ppnew.1d3b <- FeaturePlot(GEX.seur,features = c("Spp1_up.P28only"), ncol = 1,
                        #cols = c("lightgrey","red"), 
                        #cols = rev(mapal),
                        raster = T, pt.size = 3.5,
            min.cutoff = -0.1, max.cutoff = 0.54) & scale_color_gradientn(colors = rev(mapal), limits= c(-0.1,0.64), 
                                                                          breaks=c(0.0,0.2,0.4,0.6))
ppnew.1d3b

ppnew.2

DAM signature expression as feature plot (try top 50, top100, top250, top500)

DAM.sig <- read.csv("figures1002/new/ranking_of_DAM_indicator_genes.csv")
DAM.sig[1:8,]
##   锘縍anking     Gene Frequence Total.score
## 1          1     Spp1        11    35.98221
## 2          2     Apoe        11    31.44704
## 3          3   Ifitm3         9    20.82901
## 4          4 Ifi27l2a         9    20.44047
## 5          5     Lyz2        11    20.22463
## 6          6     Cst7         9    19.48372
## 7          7     Cd74         7    18.24720
## 8          8   Lgals3         7    18.24180
DAM.list <- list(top50=DAM.sig$Gene[1:50],
                 top100=DAM.sig$Gene[1:100],
                 top250=DAM.sig$Gene[1:250],
                 top500=DAM.sig$Gene[1:500])
DAM.list
## $top50
##  [1] "Spp1"     "Apoe"     "Ifitm3"   "Ifi27l2a" "Lyz2"     "Cst7"    
##  [7] "Cd74"     "Lgals3"   "Lpl"      "Cd52"     "Ccl5"     "Ccl12"   
## [13] "Lgals3bp" "Cd63"     "H2-Ab1"   "Fn1"      "H2-Aa"    "Fxyd5"   
## [19] "Ccl4"     "Cybb"     "Bst2"     "Cstb"     "H2-Eb1"   "Fth1"    
## [25] "Vim"      "Tspo"     "Ctsb"     "Ccl3"     "Axl"      "Anxa5"   
## [31] "Isg15"    "Lgals1"   "Ccl2"     "Ifi204"   "Igf1"     "Ch25h"   
## [37] "Mif"      "Cxcl10"   "Plin2"    "Fabp5"    "Il1b"     "Crip1"   
## [43] "Slfn2"    "B2m"      "Irf7"     "Cd72"     "Capg"     "Ms4a6c"  
## [49] "Ifit3"    "Apoc1"   
## 
## $top100
##   [1] "Spp1"     "Apoe"     "Ifitm3"   "Ifi27l2a" "Lyz2"     "Cst7"    
##   [7] "Cd74"     "Lgals3"   "Lpl"      "Cd52"     "Ccl5"     "Ccl12"   
##  [13] "Lgals3bp" "Cd63"     "H2-Ab1"   "Fn1"      "H2-Aa"    "Fxyd5"   
##  [19] "Ccl4"     "Cybb"     "Bst2"     "Cstb"     "H2-Eb1"   "Fth1"    
##  [25] "Vim"      "Tspo"     "Ctsb"     "Ccl3"     "Axl"      "Anxa5"   
##  [31] "Isg15"    "Lgals1"   "Ccl2"     "Ifi204"   "Igf1"     "Ch25h"   
##  [37] "Mif"      "Cxcl10"   "Plin2"    "Fabp5"    "Il1b"     "Crip1"   
##  [43] "Slfn2"    "B2m"      "Irf7"     "Cd72"     "Capg"     "Ms4a6c"  
##  [49] "Ifit3"    "Apoc1"    "Fcgr4"    "Il1rn"    "Wfdc17"   "Cxcl2"   
##  [55] "Cxcl16"   "Prdx1"    "Rtp4"     "H2-D1"    "Pkm"      "Stat1"   
##  [61] "Anxa2"    "Gpnmb"    "Zbp1"     "Ftl1"     "Ldha"     "Npc2"    
##  [67] "Ly6a"     "Oasl2"    "Gapdh"    "Prdx5"    "Gbp2"     "Grn"     
##  [73] "Ifi207"   "Ifitm2"   "Tlr2"     "Txn1"     "Phf11b"   "Ctsz"    
##  [79] "Pianp"    "Cd36"     "Itgax"    "Fgl2"     "Ccl6"     "Iigp1"   
##  [85] "Ccl7"     "H2-K1"    "Pld3"     "Tpi1"     "Akr1a1"   "Usp18"   
##  [91] "Rab7b"    "Tmsb10"   "Cd44"     "Aldoa"    "Hcar2"    "Acod1"   
##  [97] "Cd300lb"  "Ctsc"     "Gpr65"    "H2-Q7"   
## 
## $top250
##   [1] "Spp1"     "Apoe"     "Ifitm3"   "Ifi27l2a" "Lyz2"     "Cst7"    
##   [7] "Cd74"     "Lgals3"   "Lpl"      "Cd52"     "Ccl5"     "Ccl12"   
##  [13] "Lgals3bp" "Cd63"     "H2-Ab1"   "Fn1"      "H2-Aa"    "Fxyd5"   
##  [19] "Ccl4"     "Cybb"     "Bst2"     "Cstb"     "H2-Eb1"   "Fth1"    
##  [25] "Vim"      "Tspo"     "Ctsb"     "Ccl3"     "Axl"      "Anxa5"   
##  [31] "Isg15"    "Lgals1"   "Ccl2"     "Ifi204"   "Igf1"     "Ch25h"   
##  [37] "Mif"      "Cxcl10"   "Plin2"    "Fabp5"    "Il1b"     "Crip1"   
##  [43] "Slfn2"    "B2m"      "Irf7"     "Cd72"     "Capg"     "Ms4a6c"  
##  [49] "Ifit3"    "Apoc1"    "Fcgr4"    "Il1rn"    "Wfdc17"   "Cxcl2"   
##  [55] "Cxcl16"   "Prdx1"    "Rtp4"     "H2-D1"    "Pkm"      "Stat1"   
##  [61] "Anxa2"    "Gpnmb"    "Zbp1"     "Ftl1"     "Ldha"     "Npc2"    
##  [67] "Ly6a"     "Oasl2"    "Gapdh"    "Prdx5"    "Gbp2"     "Grn"     
##  [73] "Ifi207"   "Ifitm2"   "Tlr2"     "Txn1"     "Phf11b"   "Ctsz"    
##  [79] "Pianp"    "Cd36"     "Itgax"    "Fgl2"     "Ccl6"     "Iigp1"   
##  [85] "Ccl7"     "H2-K1"    "Pld3"     "Tpi1"     "Akr1a1"   "Usp18"   
##  [91] "Rab7b"    "Tmsb10"   "Cd44"     "Aldoa"    "Hcar2"    "Acod1"   
##  [97] "Cd300lb"  "Ctsc"     "Gpr65"    "H2-Q7"    "Cdkn1a"   "Psat1"   
## [103] "Trim30a"  "Cxcl14"   "Srgn"     "Mmp12"    "Bcl2a1b"  "Tap1"    
## [109] "AB124611" "Xaf1"     "Ly6e"     "Psme1"    "Ctsl"     "Sp100"   
## [115] "Cxcr4"    "Psmb8"    "AA467197" "Pgk1"     "Emp3"     "Csf1"    
## [121] "Mcemp1"   "Gusb"     "Pim1"     "Ctse"     "Cox4i1"   "Il12b"   
## [127] "Msrb1"    "Npm1"     "Alcam"    "Psme2"    "Apoc2"    "Bhlhe40" 
## [133] "Stmn1"    "Gm4951"   "Pla2g7"   "Plaur"    "Tor3a"    "Hspe1"   
## [139] "Tpt1"     "Ndufa1"   "Flt1"     "Tgfbi"    "Cox6c"    "Irgm1"   
## [145] "Ifit2"    "Uba52"    "Igf2r"    "Bola2"    "Ank"      "Tyrobp"  
## [151] "Tpm4"     "Ass1"     "Ms4a4c"   "Ifit1"    "Ybx1"     "Sod2"    
## [157] "Cox8a"    "Fam20c"   "Oas1a"    "Arg1"     "Ms4a7"    "Smpdl3a" 
## [163] "Sh3pxd2b" "Fau"      "Gnas"     "Phf11d"   "Ehd1"     "Saa3"    
## [169] "Cox5a"    "Atox1"    "Id3"      "Lrpap1"   "Olr1"     "Sh3bgrl3"
## [175] "Sash1"    "Hint1"    "Il2rg"    "Ctsd"     "Postn"    "H2-T23"  
## [181] "Ucp2"     "Sdc3"     "Uqcrq"    "Cox6a1"   "Cd300lf"  "Syngr1"  
## [187] "Timp2"    "Atp5e"    "Id2"      "S100a6"   "Cd14"     "Tubb6"   
## [193] "Anp32b"   "Fcgr2b"   "Cd83"     "Psmb9"    "Bcl2a1a"  "Aprt"    
## [199] "Mfsd12"   "Psap"     "Cox6b1"   "Lilr4b"   "Plac8"    "Glrx"    
## [205] "Scimp"    "Rilpl2"   "Psmb6"    "Gm11808"  "Chmp4b"   "Actr3b"  
## [211] "Ly86"     "Fundc2"   "Ifi211"   "Hif1a"    "Serf2"    "Dram1"   
## [217] "Parp14"   "Ptgs2"    "Cxcl9"    "Myo5a"    "Gabarap"  "Cyp4f18" 
## [223] "Shisa5"   "Lilrb4a"  "Cpd"      "Iqgap1"   "Slfn5"    "Tnfaip2" 
## [229] "Nme1"     "Cotl1"    "Tagln2"   "Mt1"      "Mpeg1"    "C3"      
## [235] "Ube2l6"   "Clic4"    "Naaa"     "Gas7"     "Sdcbp"    "Nampt"   
## [241] "Ell2"     "Samhd1"   "Rtcb"     "Eef1g"    "Hmgb2"    "Gng5"    
## [247] "Nfil3"    "Adora1"   "Tpd52"    "Ptger4"  
## 
## $top500
##   [1] "Spp1"          "Apoe"          "Ifitm3"        "Ifi27l2a"     
##   [5] "Lyz2"          "Cst7"          "Cd74"          "Lgals3"       
##   [9] "Lpl"           "Cd52"          "Ccl5"          "Ccl12"        
##  [13] "Lgals3bp"      "Cd63"          "H2-Ab1"        "Fn1"          
##  [17] "H2-Aa"         "Fxyd5"         "Ccl4"          "Cybb"         
##  [21] "Bst2"          "Cstb"          "H2-Eb1"        "Fth1"         
##  [25] "Vim"           "Tspo"          "Ctsb"          "Ccl3"         
##  [29] "Axl"           "Anxa5"         "Isg15"         "Lgals1"       
##  [33] "Ccl2"          "Ifi204"        "Igf1"          "Ch25h"        
##  [37] "Mif"           "Cxcl10"        "Plin2"         "Fabp5"        
##  [41] "Il1b"          "Crip1"         "Slfn2"         "B2m"          
##  [45] "Irf7"          "Cd72"          "Capg"          "Ms4a6c"       
##  [49] "Ifit3"         "Apoc1"         "Fcgr4"         "Il1rn"        
##  [53] "Wfdc17"        "Cxcl2"         "Cxcl16"        "Prdx1"        
##  [57] "Rtp4"          "H2-D1"         "Pkm"           "Stat1"        
##  [61] "Anxa2"         "Gpnmb"         "Zbp1"          "Ftl1"         
##  [65] "Ldha"          "Npc2"          "Ly6a"          "Oasl2"        
##  [69] "Gapdh"         "Prdx5"         "Gbp2"          "Grn"          
##  [73] "Ifi207"        "Ifitm2"        "Tlr2"          "Txn1"         
##  [77] "Phf11b"        "Ctsz"          "Pianp"         "Cd36"         
##  [81] "Itgax"         "Fgl2"          "Ccl6"          "Iigp1"        
##  [85] "Ccl7"          "H2-K1"         "Pld3"          "Tpi1"         
##  [89] "Akr1a1"        "Usp18"         "Rab7b"         "Tmsb10"       
##  [93] "Cd44"          "Aldoa"         "Hcar2"         "Acod1"        
##  [97] "Cd300lb"       "Ctsc"          "Gpr65"         "H2-Q7"        
## [101] "Cdkn1a"        "Psat1"         "Trim30a"       "Cxcl14"       
## [105] "Srgn"          "Mmp12"         "Bcl2a1b"       "Tap1"         
## [109] "AB124611"      "Xaf1"          "Ly6e"          "Psme1"        
## [113] "Ctsl"          "Sp100"         "Cxcr4"         "Psmb8"        
## [117] "AA467197"      "Pgk1"          "Emp3"          "Csf1"         
## [121] "Mcemp1"        "Gusb"          "Pim1"          "Ctse"         
## [125] "Cox4i1"        "Il12b"         "Msrb1"         "Npm1"         
## [129] "Alcam"         "Psme2"         "Apoc2"         "Bhlhe40"      
## [133] "Stmn1"         "Gm4951"        "Pla2g7"        "Plaur"        
## [137] "Tor3a"         "Hspe1"         "Tpt1"          "Ndufa1"       
## [141] "Flt1"          "Tgfbi"         "Cox6c"         "Irgm1"        
## [145] "Ifit2"         "Uba52"         "Igf2r"         "Bola2"        
## [149] "Ank"           "Tyrobp"        "Tpm4"          "Ass1"         
## [153] "Ms4a4c"        "Ifit1"         "Ybx1"          "Sod2"         
## [157] "Cox8a"         "Fam20c"        "Oas1a"         "Arg1"         
## [161] "Ms4a7"         "Smpdl3a"       "Sh3pxd2b"      "Fau"          
## [165] "Gnas"          "Phf11d"        "Ehd1"          "Saa3"         
## [169] "Cox5a"         "Atox1"         "Id3"           "Lrpap1"       
## [173] "Olr1"          "Sh3bgrl3"      "Sash1"         "Hint1"        
## [177] "Il2rg"         "Ctsd"          "Postn"         "H2-T23"       
## [181] "Ucp2"          "Sdc3"          "Uqcrq"         "Cox6a1"       
## [185] "Cd300lf"       "Syngr1"        "Timp2"         "Atp5e"        
## [189] "Id2"           "S100a6"        "Cd14"          "Tubb6"        
## [193] "Anp32b"        "Fcgr2b"        "Cd83"          "Psmb9"        
## [197] "Bcl2a1a"       "Aprt"          "Mfsd12"        "Psap"         
## [201] "Cox6b1"        "Lilr4b"        "Plac8"         "Glrx"         
## [205] "Scimp"         "Rilpl2"        "Psmb6"         "Gm11808"      
## [209] "Chmp4b"        "Actr3b"        "Ly86"          "Fundc2"       
## [213] "Ifi211"        "Hif1a"         "Serf2"         "Dram1"        
## [217] "Parp14"        "Ptgs2"         "Cxcl9"         "Myo5a"        
## [221] "Gabarap"       "Cyp4f18"       "Shisa5"        "Lilrb4a"      
## [225] "Cpd"           "Iqgap1"        "Slfn5"         "Tnfaip2"      
## [229] "Nme1"          "Cotl1"         "Tagln2"        "Mt1"          
## [233] "Mpeg1"         "C3"            "Ube2l6"        "Clic4"        
## [237] "Naaa"          "Gas7"          "Sdcbp"         "Nampt"        
## [241] "Ell2"          "Samhd1"        "Rtcb"          "Eef1g"        
## [245] "Hmgb2"         "Gng5"          "Nfil3"         "Adora1"       
## [249] "Tpd52"         "Ptger4"        "Eif2ak2"       "Areg"         
## [253] "Rsad2"         "Slc31a1"       "Gm2000"        "Cox7c"        
## [257] "Tmem256"       "Cox5b"         "Cyba"          "Il18bp"       
## [261] "Selenow"       "Myl6"          "H2-DMa"        "Rai14"        
## [265] "Dab2"          "Hmox1"         "Tmsb4x"        "Ndufa2"       
## [269] "Cd68"          "Ccdc86"        "Atp6v1a"       "Cox7a2"       
## [273] "Calm1"         "Uqcrh"         "Socs2"         "Gpi1"         
## [277] "Ubl5"          "Colec12"       "Ifit3b"        "Scpep1"       
## [281] "S100a11"       "F3"            "Ms4a6d"        "Hacd2"        
## [285] "Chil3"         "Tuba1c"        "Rap2b"         "Gp49a"        
## [289] "Milr1"         "Fcgr1"         "Stat2"         "Atp5j2"       
## [293] "Uqcr10"        "Ssr4"          "Bcar3"         "Gm49368"      
## [297] "Tmem106a"      "Tubb5"         "Naca"          "Hspa8"        
## [301] "Atp5k"         "Bag1"          "Sec61b"        "Gyg"          
## [305] "Cox7b"         "Ly6c2"         "Top2a"         "Aldh2"        
## [309] "Dpp7"          "Eif3k"         "Cd48"          "C4b"          
## [313] "Pycard"        "Atp5l"         "Uqcr11"        "Osm"          
## [317] "Prelid1"       "Rnf213"        "Prdx4"         "Arpc1b"       
## [321] "Ndufv3"        "Sp110"         "Ndufa13"       "Abca1"        
## [325] "Gm1673"        "Wdr89"         "Sh3bgrl"       "Smdt1"        
## [329] "Gatm"          "Gpr84"         "Slamf8"        "Ccrl2"        
## [333] "Tomm7"         "Gas8"          "Ly6i"          "Bcl2a1d"      
## [337] "Esd"           "Dhx58"         "Ctsa"          "Rxrg"         
## [341] "Prdx6"         "Ndufc1"        "Polr2f"        "Sdc4"         
## [345] "Atp5j"         "Litaf"         "Atp6v0e"       "Cspg4"        
## [349] "Ranbp1"        "Ifi35"         "Fcer1g"        "Calm3"        
## [353] "Rhoc"          "Pde4b"         "Atp5g1"        "Gpx3"         
## [357] "Psmb1"         "Gpx1"          "Eef1a1"        "Ly9"          
## [361] "Igtp"          "H2-Q6"         "Herc6"         "Cd9"          
## [365] "Ran"           "Hebp1"         "Eno1"          "Cdk18"        
## [369] "Eef1b2"        "Gbp7"          "Glipr1"        "Atp6v1f"      
## [373] "H2-DMb1"       "Btf3"          "Slc25a3"       "Brdt"         
## [377] "Csf2ra"        "Eif3f"         "Hpse"          "Gm14305"      
## [381] "Gm14410"       "H2-Q4"         "Ndufb9"        "Epsti1"       
## [385] "Dnaaf3"        "Pf4"           "S100a4"        "Atp5h"        
## [389] "Apobec3"       "Hsp90ab1"      "Tnf"           "Ctss"         
## [393] "Gas6"          "Gbp3"          "Gng12"         "Nceh1"        
## [397] "Mgst1"         "Cfl1"          "Dtnbp1"        "Slc2a6"       
## [401] "Eif5a"         "Atp5c1"        "ptchd1"        "Ptchd1"       
## [405] "Psma7"         "Lap3"          "Rbm3"          "Cycs"         
## [409] "Cox6a2"        "Abcg1"         "Prr15"         "Ahnak"        
## [413] "Ndufb7"        "Nrp1"          "Lmnb1"         "Ninj1"        
## [417] "Mrpl33"        "Arpc3"         "Phlda1"        "Acp5"         
## [421] "Atp5g3"        "C5ar1"         "Arl5c"         "Parp9"        
## [425] "Ndufb8"        "Txndc17"       "Tbca"          "Gm49339"      
## [429] "Tma7"          "Fblim1"        "Eif3h"         "Psme2b"       
## [433] "Mrpl30"        "Arpc2"         "Ptma"          "Rac2"         
## [437] "Ctsh"          "Dcstamp"       "Clec4n"        "Dbi"          
## [441] "H2-T22"        "Sem1"          "Msr1"          "Bax"          
## [445] "S100a10"       "Fabp3"         "Snrpg"         "Bcl2a1"       
## [449] "Serpine2"      "Snrpd2"        "Cd180"         "Pgam1"        
## [453] "2310061I04Rik" "S100a1"        "Serpinb6a"     "Actg1"        
## [457] "Uqcc2"         "Tuba1b"        "Neurl2"        "Gcnt2"        
## [461] "Smc2"          "Atp5b"         "Ifih1"         "Nap1l1"       
## [465] "Cd274"         "Rgs1"          "Parp12"        "Serpine1"     
## [469] "Nsa2"          "Cebpb"         "Csf2rb"        "Cfb"          
## [473] "Ndufa6"        "Sdf2l1"        "Anxa4"         "Psma5"        
## [477] "Nfkbiz"        "Ctla2a"        "Atp5o"         "Ube2l3"       
## [481] "Lipa"          "Pfn1"          "Ndufa7"        "Ndufs6"       
## [485] "Psmb10"        "Mapkapk2"      "Aif1"          "Smc4"         
## [489] "Itgb1bp1"      "Nme2"          "Pfdn5"         "Rassf3"       
## [493] "1810058I24Rik" "Pgls"          "Clec4d"        "Mrpl54"       
## [497] "Tmem160"       "Pomp"          "Slc7a11"       "F10"
lapply(DAM.list, length)
## $top50
## [1] 50
## 
## $top100
## [1] 100
## 
## $top250
## [1] 250
## 
## $top500
## [1] 500
GEX.seur <- add_geneset_score(GEX.seur, DAM.list$top50, "DAM.sig_top50")
## Summarizing data
GEX.seur <- add_geneset_score(GEX.seur, DAM.list$top100, "DAM.sig_top100")
## Summarizing data
GEX.seur <- add_geneset_score(GEX.seur, DAM.list$top250, "DAM.sig_top250")
## Summarizing data
GEX.seur <- add_geneset_score(GEX.seur, DAM.list$top500, "DAM.sig_top500")
## Summarizing data
ppnew.2a <- FeaturePlot(GEX.seur,features = c("DAM.sig_top50","DAM.sig_top100","DAM.sig_top250","DAM.sig_top500"), ncol = 4,
                        raster = T, pt.size = 3.5,
            keep.scale = "all")
ppnew.2a

ppnew.2b <- FeaturePlot(GEX.seur,features = c("DAM.sig_top50","DAM.sig_top100","DAM.sig_top250","DAM.sig_top500"), ncol = 4,
                        raster = T, pt.size = 3.5,
                        cols = c("lightgrey","red"),
            keep.scale = "all")
ppnew.2b

ppnew.2c <- FeaturePlot(GEX.seur,features = c("DAM.sig_top50","DAM.sig_top100","DAM.sig_top250","DAM.sig_top500"), ncol = 4,
                        raster = T, pt.size = 3.5,
            keep.scale = "all") & scale_color_gradientn(colors = rev(mapal), limits= c(-0.4,1.65), breaks = c(0.0,0.5,1.0,1.5))
ppnew.2c

ppnew.2d <- FeaturePlot(GEX.seur,features = c("DAM.sig_top50","DAM.sig_top100","DAM.sig_top250","DAM.sig_top500"), ncol = 4,
                        raster = T, pt.size = 3.5,
            keep.scale = "all") & scale_color_gradientn(colors = rev(mapal))
ppnew.2d

comparison

each cluster, compare CTRvsGFP

head(GEX.seur@meta.data)
##                     orig.ident nCount_RNA nFeature_RNA percent.mt percent.rb
## AAACCCAAGAAATCCA-1 plus_SIMPLE       2998         1458   2.701801   15.01001
## AAACCCAAGCGATCGA-1 plus_SIMPLE       5474         2107   5.608330   17.97589
## AAACCCAAGGTAGCCA-1 plus_SIMPLE       7120         2540   3.609551   12.48596
## AAACCCAAGGTTGAGC-1 plus_SIMPLE       4650         1705   4.430108   19.03226
## AAACCCAAGTCACTCA-1 plus_SIMPLE       2150          894   2.790698   10.04651
## AAACCCAAGTGGAATT-1 plus_SIMPLE       3738         1568   7.276619   20.43874
##                      FB.info     S.Score   G2M.Score Phase RNA_snn_res.1.5
## AAACCCAAGAAATCCA-1 P05.GFP.1 -0.06682298 -0.09221214    G1               2
## AAACCCAAGCGATCGA-1 P05.MIG.2 -0.09983915 -0.03784219    G1               1
## AAACCCAAGGTAGCCA-1 P28.GFP.1 -0.07149868 -0.02169489    G1               9
## AAACCCAAGGTTGAGC-1 P28.GFP.1  0.03069876 -0.11000513     S               5
## AAACCCAAGTCACTCA-1 P28.GFP.1 -0.01015567 -0.02782644    G1              11
## AAACCCAAGTGGAATT-1 P05.MIG.2 -0.02148368 -0.04812337    G1               6
##                    seurat_clusters RNA_snn_res.1.2     cnt sort_clusters
## AAACCCAAGAAATCCA-1               1               1 P05.GFP             1
## AAACCCAAGCGATCGA-1               1               1 P05.MIG             1
## AAACCCAAGGTAGCCA-1               5               5 P28.GFP             5
## AAACCCAAGGTTGAGC-1               4               4 P28.GFP             4
## AAACCCAAGTCACTCA-1               7               7 P28.GFP             7
## AAACCCAAGTGGAATT-1              10              10 P05.MIG            10
##                    pANN_0.25_0.005_1267 DoubletFinder0.05 pANN_0.25_0.005_2534
## AAACCCAAGAAATCCA-1           0.00000000           Singlet           0.00000000
## AAACCCAAGCGATCGA-1           0.49704142           Singlet           0.49112426
## AAACCCAAGGTAGCCA-1           0.44378698           Singlet           0.42603550
## AAACCCAAGGTTGAGC-1           0.14201183           Singlet           0.11834320
## AAACCCAAGTCACTCA-1           0.00000000           Singlet           0.00000000
## AAACCCAAGTGGAATT-1           0.01183432           Singlet           0.01775148
##                    DoubletFinder0.1 preAnno1 preAnno2 preAnno    FB.new  Anno1
## AAACCCAAGAAATCCA-1          Singlet   MIG.a2    MIG.a     MIG P05.GFP.1 MIG.a2
## AAACCCAAGCGATCGA-1          Doublet   MIG.a2    MIG.a     MIG P05.MIG.2 MIG.a2
## AAACCCAAGGTAGCCA-1          Doublet   MIG.b7    MIG.b     MIG P28.GFP.1 MIG.b7
## AAACCCAAGGTTGAGC-1          Singlet   MIG.b7    MIG.b     MIG P28.GFP.1 MIG.b3
## AAACCCAAGTCACTCA-1          Singlet   MIG.b6    MIG.b     MIG P28.GFP.1 MIG.b4
## AAACCCAAGTGGAATT-1          Singlet   MIG.a3    MIG.a     MIG P05.MIG.2 MIG.a4
##                    Anno2  rep Spp1_up.P06only Spp1_up.both Spp1_up.P28only
## AAACCCAAGAAATCCA-1 MIG.a rep1      0.07799678   0.01572847       0.2361100
## AAACCCAAGCGATCGA-1 MIG.a rep2     -0.04681886  -0.05326058       0.2947471
## AAACCCAAGGTAGCCA-1 MIG.b rep1     -0.03821920   0.29210607       0.4436883
## AAACCCAAGGTTGAGC-1 MIG.b rep1     -0.00965059   0.10195132       0.4066607
## AAACCCAAGTCACTCA-1 MIG.b rep1      0.03137777   0.29949531       0.3521602
## AAACCCAAGTGGAATT-1 MIG.a rep2     -0.02037658  -0.01007519       0.3301183
##                    DAM.sig_top50 DAM.sig_top100 DAM.sig_top250 DAM.sig_top500
## AAACCCAAGAAATCCA-1   -0.15072706    -0.18592966    -0.13817389     0.02798874
## AAACCCAAGCGATCGA-1   -0.23592021    -0.16868345    -0.06001030     0.12662573
## AAACCCAAGGTAGCCA-1    0.79968333     0.61875031     0.48767439     0.50833408
## AAACCCAAGGTTGAGC-1    0.25948931     0.20225591     0.17024206     0.31339640
## AAACCCAAGTCACTCA-1    0.90775071     0.64612166     0.38459858     0.38874980
## AAACCCAAGTGGAATT-1   -0.02039969    -0.04251335    -0.02195155     0.17200874
##                    Cycling_score
## AAACCCAAGAAATCCA-1  -0.072838959
## AAACCCAAGCGATCGA-1  -0.045974881
## AAACCCAAGGTAGCCA-1  -0.006387924
## AAACCCAAGGTTGAGC-1  -0.010882599
## AAACCCAAGTCACTCA-1  -0.036946642
## AAACCCAAGTGGAATT-1  -0.007985408
ppnew.2.v1 <- VlnPlot(GEX.seur, features = c("DAM.sig_top50","DAM.sig_top100","DAM.sig_top250","DAM.sig_top500"), 
                      ncol = 2, pt.size = 0, raster = F) & geom_jitter(alpha=0.15, shape=16, width = 0.2, size = 0.02)
ppnew.2.v1

ppnew.2.v1 <- VlnPlot(GEX.seur, features = c("DAM.sig_top50","DAM.sig_top100","DAM.sig_top250","DAM.sig_top500"), 
                      same.y.lims = T,
                      ncol = 2, pt.size = 0, raster = F) & geom_jitter(alpha=0.15, shape=16, width = 0.2, size = 0.02)
ppnew.2.v1

ppnew.2.v1 <- VlnPlot(GEX.seur, features = c("DAM.sig_top50","DAM.sig_top100","DAM.sig_top250","DAM.sig_top500"), 
                      group.by = "cnt", cols = color.cnt,
                      ncol = 2, pt.size = 0, raster = F) & geom_jitter(alpha=0.15, shape=16, width = 0.2, size = 0.02)
ppnew.2.v1

ppnew.2.v1 <- VlnPlot(GEX.seur, features = c("DAM.sig_top50","DAM.sig_top100","DAM.sig_top250","DAM.sig_top500"), 
                      same.y.lims = T,
                      group.by = "cnt", cols = color.cnt,
                      ncol = 2, pt.size = 0, raster = F) & 
    geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
  stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) &
  ggpubr::stat_compare_means(aes(lable = ..p.signif..), 
                               method = "wilcox.test",
                               comparisons = list(c("P05.CTR","P05.GFP"),
                                                  c("P28.CTR","P28.GFP")),
                               label.y = c(1,1.4),
                               size=3
                               )
ppnew.2.v1

ppnew.2.v2 <- VlnPlot(GEX.seur, features = c("DAM.sig_top50","DAM.sig_top100","DAM.sig_top250","DAM.sig_top500"), 
                      split.by = "cnt", cols = color.cnt,
                      ncol = 1, pt.size = 0, raster = F) #& geom_jitter(alpha=0.15, shape=16, width = 0.2, size = 0.02)
ppnew.2.v2

subset

subset ‘MIG.a - P05’ and ‘MIG.b - P28’ to check ‘CTR vs GFP’ in each cluster

GEX.P05_MIGa  <- subset(GEX.seur, subset = Anno2 == "MIG.a" & cnt %in% c("P05.CTR","P05.MIG","P05.GFP"))
GEX.P05_MIGa
## An object of class Seurat 
## 18371 features across 11960 samples within 1 assay 
## Active assay: RNA (18371 features, 1389 variable features)
##  3 dimensional reductions calculated: pca, tsne, umap
GEX.P28_MIGb  <- subset(GEX.seur, subset = Anno2 == "MIG.b" & cnt %in% c("P28.CTR","P28.MIG","P28.GFP"))
GEX.P28_MIGb
## An object of class Seurat 
## 18371 features across 11099 samples within 1 assay 
## Active assay: RNA (18371 features, 1389 variable features)
##  3 dimensional reductions calculated: pca, tsne, umap
cnta
sort(as.character(unique(GEX.P05_MIGa$Anno1)))
## [1] "MIG.a1" "MIG.a2" "MIG.a3" "MIG.a4" "MIG.a5" "MIG.a6" "MIG.a7" "MIG.a8"
sort(as.character(unique(GEX.P05_MIGa$cnt)))
## [1] "P05.CTR" "P05.GFP" "P05.MIG"
cnta <- sapply(sort(as.character(unique(GEX.P05_MIGa$Anno1))),function(x){
  paste0(x,"_",sort(as.character(unique(GEX.P05_MIGa$cnt)))[c(1,3,2)])
})
cnta <- as.vector(unlist(cnta))
cnta
##  [1] "MIG.a1_P05.CTR" "MIG.a1_P05.MIG" "MIG.a1_P05.GFP" "MIG.a2_P05.CTR"
##  [5] "MIG.a2_P05.MIG" "MIG.a2_P05.GFP" "MIG.a3_P05.CTR" "MIG.a3_P05.MIG"
##  [9] "MIG.a3_P05.GFP" "MIG.a4_P05.CTR" "MIG.a4_P05.MIG" "MIG.a4_P05.GFP"
## [13] "MIG.a5_P05.CTR" "MIG.a5_P05.MIG" "MIG.a5_P05.GFP" "MIG.a6_P05.CTR"
## [17] "MIG.a6_P05.MIG" "MIG.a6_P05.GFP" "MIG.a7_P05.CTR" "MIG.a7_P05.MIG"
## [21] "MIG.a7_P05.GFP" "MIG.a8_P05.CTR" "MIG.a8_P05.MIG" "MIG.a8_P05.GFP"
ctta <- list()
for(i in 1:8){
  ctta[[i]] <- cnta[c(i*3-2,i*3)]
}
ctta
## [[1]]
## [1] "MIG.a1_P05.CTR" "MIG.a1_P05.GFP"
## 
## [[2]]
## [1] "MIG.a2_P05.CTR" "MIG.a2_P05.GFP"
## 
## [[3]]
## [1] "MIG.a3_P05.CTR" "MIG.a3_P05.GFP"
## 
## [[4]]
## [1] "MIG.a4_P05.CTR" "MIG.a4_P05.GFP"
## 
## [[5]]
## [1] "MIG.a5_P05.CTR" "MIG.a5_P05.GFP"
## 
## [[6]]
## [1] "MIG.a6_P05.CTR" "MIG.a6_P05.GFP"
## 
## [[7]]
## [1] "MIG.a7_P05.CTR" "MIG.a7_P05.GFP"
## 
## [[8]]
## [1] "MIG.a8_P05.CTR" "MIG.a8_P05.GFP"
GEX.P05_MIGa$Anno1_cnt <- factor(paste0(as.character(GEX.P05_MIGa$Anno1),
                                        "_",
                                        as.character(GEX.P05_MIGa$cnt)),
                                 levels = cnta)
cntb
sort(as.character(unique(GEX.P28_MIGb$Anno1)))
## [1] "MIG.b1" "MIG.b2" "MIG.b3" "MIG.b4" "MIG.b5" "MIG.b6" "MIG.b7" "MIG.b8"
sort(as.character(unique(GEX.P28_MIGb$cnt)))
## [1] "P28.CTR" "P28.GFP" "P28.MIG"
cntb <- sapply(sort(as.character(unique(GEX.P28_MIGb$Anno1))),function(x){
  paste0(x,"_",sort(as.character(unique(GEX.P28_MIGb$cnt)))[c(1,3,2)])
})
cntb <- as.vector(unlist(cntb))
cntb
##  [1] "MIG.b1_P28.CTR" "MIG.b1_P28.MIG" "MIG.b1_P28.GFP" "MIG.b2_P28.CTR"
##  [5] "MIG.b2_P28.MIG" "MIG.b2_P28.GFP" "MIG.b3_P28.CTR" "MIG.b3_P28.MIG"
##  [9] "MIG.b3_P28.GFP" "MIG.b4_P28.CTR" "MIG.b4_P28.MIG" "MIG.b4_P28.GFP"
## [13] "MIG.b5_P28.CTR" "MIG.b5_P28.MIG" "MIG.b5_P28.GFP" "MIG.b6_P28.CTR"
## [17] "MIG.b6_P28.MIG" "MIG.b6_P28.GFP" "MIG.b7_P28.CTR" "MIG.b7_P28.MIG"
## [21] "MIG.b7_P28.GFP" "MIG.b8_P28.CTR" "MIG.b8_P28.MIG" "MIG.b8_P28.GFP"
cttb <- list()
for(i in 1:8){
  cttb[[i]] <- cntb[c(i*3-2,i*3)]
}
cttb
## [[1]]
## [1] "MIG.b1_P28.CTR" "MIG.b1_P28.GFP"
## 
## [[2]]
## [1] "MIG.b2_P28.CTR" "MIG.b2_P28.GFP"
## 
## [[3]]
## [1] "MIG.b3_P28.CTR" "MIG.b3_P28.GFP"
## 
## [[4]]
## [1] "MIG.b4_P28.CTR" "MIG.b4_P28.GFP"
## 
## [[5]]
## [1] "MIG.b5_P28.CTR" "MIG.b5_P28.GFP"
## 
## [[6]]
## [1] "MIG.b6_P28.CTR" "MIG.b6_P28.GFP"
## 
## [[7]]
## [1] "MIG.b7_P28.CTR" "MIG.b7_P28.GFP"
## 
## [[8]]
## [1] "MIG.b8_P28.CTR" "MIG.b8_P28.GFP"
GEX.P28_MIGb$Anno1_cnt <- factor(paste0(as.character(GEX.P28_MIGb$Anno1),
                                        "_",
                                        as.character(GEX.P28_MIGb$cnt)),
                                 levels = cntb)
newplot
ppnew.2.v2a <- VlnPlot(GEX.P05_MIGa, features = c("DAM.sig_top50","DAM.sig_top100","DAM.sig_top250","DAM.sig_top500"), 
                      same.y.lims = T,
                      group.by = "Anno1_cnt", cols = rep(color.cnt[1:3],8),
                      ncol = 2, pt.size = 0, raster = F) & 
    geom_boxplot(outlier.size = 0, fill="white", width=0.35, size=0.15, alpha=0.55) &
  stat_summary(fun=mean, geom="point", shape=18, size=2, color="black", alpha=0.45) &
  ggpubr::stat_compare_means(aes(lable = ..p.signif..), 
                               method = "wilcox.test",
                               comparisons = ctta,
                               label.y = rep(0.53,8),
                               size=2.1
                               ) & theme(axis.text.x = element_text(size=7.5)) & ylim(c(-0.4,0.6))
ppnew.2.v2a

ppnew.2.v2b <- VlnPlot(GEX.P28_MIGb, features = c("DAM.sig_top50","DAM.sig_top100","DAM.sig_top250","DAM.sig_top500"), 
                      same.y.lims = T,
                      group.by = "Anno1_cnt", cols = c(rep(color.cnt[c(4,6,5)],7),color.cnt[6:5]),
                      ncol = 2, pt.size = 0, raster = F) & 
    geom_boxplot(outlier.size = 0, fill="white", width=0.35, size=0.15, alpha=0.55) &
  stat_summary(fun=mean, geom="point", shape=18, size=2, color="black", alpha=0.45) &
  ggpubr::stat_compare_means(aes(lable = ..p.signif..), 
                               method = "wilcox.test",
                               comparisons = cttb[c(1:5)],
                               label.y = rep(1.55,8),
                               size=2.1
                               ) & theme(axis.text.x = element_text(size=7.5)) & ylim(c(-0.5,1.75))
ppnew.2.v2b

ppnew.3

New dotplot (gene list attached)

new.marker <- as.vector(unlist(read.csv("figures1002/new/dotplot_gene.csv")))
new.marker <- unique(new.marker)
new.marker
##  [1] "Tmem119" "P2ry12"  "P2ry13"  "Cx3cr1"  "Crybb1"  "Sall1"   "Mafb"   
##  [8] "Igf1"    "Apoe"    "Cd63"    "Abcg1"   "Lpl"     "Psat1"   "Csf1"   
## [15] "Lilrb4a" "Ftl1"    "Mt1"     "Abca1"   "Gas6"    "Cxcl10"  "Ccl12"  
## [22] "Ccl2"    "Oasl2"   "Irf7"    "Oasl1"   "Tnf"     "Il1b"    "Bcl2a1b"
## [29] "Bcl2a1d" "Bcl2a1a" "Gpx4"    "Bcl2"    "Ccl5"
## mod color 
#  scale_color_gradientn(colours  = material.heat(100))

# Cell2020     
material.heat = function(n){
  colorRampPalette(
    c(
      "#283593", # indigo 800
      "#3F51B5", # indigo
      "#2196F3", # blue
      "#00BCD4", # cyan
      "#4CAF50", # green
      "#8BC34A", # light green
      "#CDDC39", # lime
      "#FFEB3B", # yellow
      "#FFC107", # amber
      "#FF9800", # orange
      "#FF5722"  # deep orange
    )
  )(n)
}

# Immunity2019, na gray
colors.Immunity <-c("#191970","#121285","#0C0C9A","#0707B0","#0101C5","#0014CF","#0033D3","#0053D8","#0072DD","#0092E1","#00B2E6",
                  "#00D1EB","#23E8CD","#7AF17B","#D2FA29","#FFEB00","#FFC300","#FF9B00","#FF8400","#FF7800","#FF6B00","#FF5F00","#FF5300",
                  "#FF4700","#F73B00","#EF2E00","#E62300","#DD1700","#D50B00","#CD0000")
ppnew.3a1 <- DotPlot(GEX.seur, features = new.marker, group.by = "Anno1"#,cols = c("midnightblue","darkorange1")
                         ) +
                 theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1)) + labs(y="Anno1")
ppnew.3a1

ppnew.3a2 <- DotPlot(GEX.seur, features = new.marker, group.by = "Anno1",
                         cols = c("midnightblue","darkorange1")) +
                 theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1)) + labs(y="Anno1")
ppnew.3a2

ppnew.3b <- DotPlot(GEX.seur, features = new.marker, group.by = "Anno1",
                         cols = c("midnightblue","darkorange1"))  +
  scale_color_gradientn(colours  = rev(mapal)) +
                 theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1)) + labs(y="Anno1")
ppnew.3b

ppnew.3c <- DotPlot(GEX.seur, features = new.marker, group.by = "Anno1",
                         cols = c("midnightblue","darkorange1")) +
  scale_color_gradientn(colours  = material.heat(100)) +
                 theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1)) + labs(y="Anno1")
ppnew.3c

ppnew.3d <- DotPlot(GEX.seur, features = new.marker, group.by = "Anno1",
                         cols = c("midnightblue","darkorange1"))  +
  scale_color_gradientn(colours  = colors.Immunity) +
                 theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1)) + labs(y="Anno1")
ppnew.3d

ppnew.4

Cycling score as feature plot

CC_gene <- Hmisc::capitalize(tolower(as.vector(unlist(cc.genes.updated.2019))))
CC_gene
##  [1] "Mcm5"     "Pcna"     "Tyms"     "Fen1"     "Mcm7"     "Mcm4"    
##  [7] "Rrm1"     "Ung"      "Gins2"    "Mcm6"     "Cdca7"    "Dtl"     
## [13] "Prim1"    "Uhrf1"    "Cenpu"    "Hells"    "Rfc2"     "Polr1b"  
## [19] "Nasp"     "Rad51ap1" "Gmnn"     "Wdr76"    "Slbp"     "Ccne2"   
## [25] "Ubr7"     "Pold3"    "Msh2"     "Atad2"    "Rad51"    "Rrm2"    
## [31] "Cdc45"    "Cdc6"     "Exo1"     "Tipin"    "Dscc1"    "Blm"     
## [37] "Casp8ap2" "Usp1"     "Clspn"    "Pola1"    "Chaf1b"   "Mrpl36"  
## [43] "E2f8"     "Hmgb2"    "Cdk1"     "Nusap1"   "Ube2c"    "Birc5"   
## [49] "Tpx2"     "Top2a"    "Ndc80"    "Cks2"     "Nuf2"     "Cks1b"   
## [55] "Mki67"    "Tmpo"     "Cenpf"    "Tacc3"    "Pimreg"   "Smc4"    
## [61] "Ccnb2"    "Ckap2l"   "Ckap2"    "Aurkb"    "Bub1"     "Kif11"   
## [67] "Anp32e"   "Tubb4b"   "Gtse1"    "Kif20b"   "Hjurp"    "Cdca3"   
## [73] "Jpt1"     "Cdc20"    "Ttk"      "Cdc25c"   "Kif2c"    "Rangap1" 
## [79] "Ncapd2"   "Dlgap5"   "Cdca2"    "Cdca8"    "Ect2"     "Kif23"   
## [85] "Hmmr"     "Aurka"    "Psrc1"    "Anln"     "Lbr"      "Ckap5"   
## [91] "Cenpe"    "Ctcf"     "Nek2"     "G2e3"     "Gas2l3"   "Cbx5"    
## [97] "Cenpa"
GEX.seur <- add_geneset_score(GEX.seur, CC_gene, "Cycling_score")
## Summarizing data
ppnew.4a <- FeaturePlot(GEX.seur,features = c("Cycling_score"), ncol = 1,
                        raster = T, pt.size = 3.5)
ppnew.4a

ppnew.4b <- FeaturePlot(GEX.seur,features = c("Cycling_score"), ncol = 1,
                        raster = T, pt.size = 3.5,
                        cols = c("lightgrey","red"))
ppnew.4b

ppnew.4c <- FeaturePlot(GEX.seur,features = c("Cycling_score"), ncol = 1,
                        raster = T, pt.size = 3.5,
                        cols = c("lightgrey","red"))& scale_color_gradientn(colors = rev(mapal), breaks = c(0.00,0.25,0.50,0.75))
ppnew.4c

#saveRDS(GEX.seur,"./GEX0811.seur.pure_Anno1.rds")

forGEO

GEX.seur <- readRDS("./GEX0811.seur.pure_Anno1.rds")
GEX.seur
## An object of class Seurat 
## 18371 features across 23692 samples within 1 assay 
## Active assay: RNA (18371 features, 1389 variable features)
##  3 dimensional reductions calculated: pca, tsne, umap
GEX.seur@meta.data[,grep("snn|pANN",colnames(GEX.seur@meta.data),value = T)] <- NULL
head(GEX.seur@meta.data)
##                     orig.ident nCount_RNA nFeature_RNA percent.mt percent.rb
## AAACCCAAGAAATCCA-1 plus_SIMPLE       2998         1458   2.701801   15.01001
## AAACCCAAGCGATCGA-1 plus_SIMPLE       5474         2107   5.608330   17.97589
## AAACCCAAGGTAGCCA-1 plus_SIMPLE       7120         2540   3.609551   12.48596
## AAACCCAAGGTTGAGC-1 plus_SIMPLE       4650         1705   4.430108   19.03226
## AAACCCAAGTCACTCA-1 plus_SIMPLE       2150          894   2.790698   10.04651
## AAACCCAAGTGGAATT-1 plus_SIMPLE       3738         1568   7.276619   20.43874
##                      FB.info     S.Score   G2M.Score Phase seurat_clusters
## AAACCCAAGAAATCCA-1 P05.GFP.1 -0.06682298 -0.09221214    G1               1
## AAACCCAAGCGATCGA-1 P05.MIG.2 -0.09983915 -0.03784219    G1               1
## AAACCCAAGGTAGCCA-1 P28.GFP.1 -0.07149868 -0.02169489    G1               5
## AAACCCAAGGTTGAGC-1 P28.GFP.1  0.03069876 -0.11000513     S               4
## AAACCCAAGTCACTCA-1 P28.GFP.1 -0.01015567 -0.02782644    G1               7
## AAACCCAAGTGGAATT-1 P05.MIG.2 -0.02148368 -0.04812337    G1              10
##                        cnt sort_clusters DoubletFinder0.05 DoubletFinder0.1
## AAACCCAAGAAATCCA-1 P05.GFP             1           Singlet          Singlet
## AAACCCAAGCGATCGA-1 P05.MIG             1           Singlet          Doublet
## AAACCCAAGGTAGCCA-1 P28.GFP             5           Singlet          Doublet
## AAACCCAAGGTTGAGC-1 P28.GFP             4           Singlet          Singlet
## AAACCCAAGTCACTCA-1 P28.GFP             7           Singlet          Singlet
## AAACCCAAGTGGAATT-1 P05.MIG            10           Singlet          Singlet
##                    preAnno1 preAnno2 preAnno    FB.new  Anno1 Anno2  rep
## AAACCCAAGAAATCCA-1   MIG.a2    MIG.a     MIG P05.GFP.1 MIG.a2 MIG.a rep1
## AAACCCAAGCGATCGA-1   MIG.a2    MIG.a     MIG P05.MIG.2 MIG.a2 MIG.a rep2
## AAACCCAAGGTAGCCA-1   MIG.b7    MIG.b     MIG P28.GFP.1 MIG.b7 MIG.b rep1
## AAACCCAAGGTTGAGC-1   MIG.b7    MIG.b     MIG P28.GFP.1 MIG.b3 MIG.b rep1
## AAACCCAAGTCACTCA-1   MIG.b6    MIG.b     MIG P28.GFP.1 MIG.b4 MIG.b rep1
## AAACCCAAGTGGAATT-1   MIG.a3    MIG.a     MIG P05.MIG.2 MIG.a4 MIG.a rep2
##                    Spp1_up.P06only Spp1_up.both Spp1_up.P28only DAM.sig_top50
## AAACCCAAGAAATCCA-1      0.07799678   0.01572847       0.2361100   -0.15072706
## AAACCCAAGCGATCGA-1     -0.04681886  -0.05326058       0.2947471   -0.23592021
## AAACCCAAGGTAGCCA-1     -0.03821920   0.29210607       0.4436883    0.79968333
## AAACCCAAGGTTGAGC-1     -0.00965059   0.10195132       0.4066607    0.25948931
## AAACCCAAGTCACTCA-1      0.03137777   0.29949531       0.3521602    0.90775071
## AAACCCAAGTGGAATT-1     -0.02037658  -0.01007519       0.3301183   -0.02039969
##                    DAM.sig_top100 DAM.sig_top250 DAM.sig_top500 Cycling_score
## AAACCCAAGAAATCCA-1    -0.18592966    -0.13817389     0.02798874  -0.072838959
## AAACCCAAGCGATCGA-1    -0.16868345    -0.06001030     0.12662573  -0.045974881
## AAACCCAAGGTAGCCA-1     0.61875031     0.48767439     0.50833408  -0.006387924
## AAACCCAAGGTTGAGC-1     0.20225591     0.17024206     0.31339640  -0.010882599
## AAACCCAAGTCACTCA-1     0.64612166     0.38459858     0.38874980  -0.036946642
## AAACCCAAGTGGAATT-1    -0.04251335    -0.02195155     0.17200874  -0.007985408
#saveRDS(GEX.seur,"I:/Shared_win/projects/202310_Spp1DAM/forGEO/seur_obj/Spp1GFP_SIMPLE.final.seur_obj.rds")